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FOREWORD 


This  report  was  prepared  by  the  Unicn  Carbide  Corporation,  Defense 
and  Space  Systems  Department,  White  Plains,  New  York,  under  Contract  DA.-U9-IU6-XZ-5II, 
Project  5710,  funded  by  the  Defense  Atomic  Support  Agency  (DASA) .  Inclusive 
dates  of  research  were  1  June  1966  to  30  June  1967.  The  report  was  submitted 
31  August  1967  by  the  AFWL  Technical  Monitor,  Captain  Guy  Spitale  (WLRP). 

The  report  has  been  divided  into  two  volumes  for  convenience .  Volume  I 
presents  the  theoretical  analysis  and  the  discussion  of  results.  Volume  II  is 
a  detailed  description  of  the  program. 

The  project  was  initiated  and  formulated  by  Dr.  C.  D.  Zerby.  The 
theoretical  derivation  was  ccmpieted  by  Dr.  H.  Brysk,  who  also  planned  and 
analyzed  the  calculations.  Programming  support  was  supplied  by  Mr.  A.  Glick 
in  writing  the  program  and  by  Mr.  E.  C.  Imperatore  in  resolving  the  systems 
problems  of  transcribing  the  tapes  containing  the  Los  Alamos  Scientific  Laboratory 
(IASL)  self -consistent-field  data  to  the  UCC  and  the  AFWL  operating  systems. 

The  program  is  written  in  FORTRAN  IV  and  is  operational  on  the  CDC  6600  computer 
at  AFWL.  We  thank  Dr.  J.  T.  Waber  (of  LASL)  who  supplied  us  the  output  of 
his  self-consistent-field  program  on  tape. 

This  report  has  been  reviewed  and  is  approved. 
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ABSTRACT 


A  computer  program  was  developed  for  the  calculation  of  photoelectric 
cross  sections,  including  angular  distributions,  using  as  input  the  results  of 
a  relativistic  Dirac-Slater  self-consistent-field  program.  The  program  was 
used  to  calculate  the  aluminum  cross  sections  over  the  range  from  1  to  150  keV 
and  uranium  cross  sections  at  four  energies  within  that  range,  and  the  results 
were  correlated  with  pre-existent  experimental  and  theoretical  data. 

Volume  I  presents  the  theoretical  analysis  and  the  discussion  of 
results.  Volume  II  is  a  detailed  description  of  the  program. 
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I.  INTRODUCTION 


This  volume  is  a  self-contained  description  of  a  program  (PELEC)  for 
the  computation  of  photoelectric  cross  sections.  The  companion  volume 
{UCC/DSSD  -  299 >  Volume  l)  develops  the  theory  used  ar.d  discusses  the  results 
obtained  with  the  program. 

Chapter-!!-  contains  the  operating  instructions.  ,  Ib-deseribes  the 

/  * 

input  cards  and  self -consistent-field  data  tape  required.  The  utilisation  of 
a  test  option  to  ascertain  the  optimum  input  parameters  in  a  series  of  runs  is 
indicated. 

Chapter  IIX  presents  sample  output.'  ■  \ 

Chapter  IV  exhibits  the  program  itself f  For  each  routine,  its  purpose 
and  the  method  of  achieving  it  are  stated.  The  routines  it  calls  and  those  it 
is  called  by  are  given,  as  well  as  the  Common  blocks  it  uses  and  the  calling 
sequence  (if  any).  The  variables  are  defined.  A  schematic  flow  chart  is 
provided.  Finally,  the  routine  is  listed. 
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II.  OPERATING  INSTRUCTIONS 


This  chapter  gives  all  the  operational  details  needed  in  order  to 
run  the  photoelectric  program.  The  input  variables  and  format  are  listed, 
along  with  the  program  diagnostics .  Tape  unit  assignments  are  given  and  the 
test  option  is  described. 

Input  Variables 

A  complete  set  of  data  consists  of  a  single  data  card.  The  variables 
described  below  a.*e  read  in  on  a  format  (9I5>  2E15.8): 

5  10  15  20  25  30  35  40  k5 _ 60 _ £>_ 

JM  KMAX  I M  NTAPE  NEDGE  IA  IB  IZ  L$0P  QV  SAVE 

Upon  completion  of  a  problem,  the  program  recycles,  reading  in  the 
next  data  card.  Tennination  of  the  program  is  obtained  by  setting  the  flag 


value  NEDGE  ■  -1. 

Name  Dimensions 

Mode 

Meaning 

JM 

I 

Maximum  order  of  Legendre  coefficient 

KMAX 

I 

Maximum  k  for  electron 

IM 

I 

Ma  dmum  l  for  photon 

NTAPE 

I 

The  logical  tape  unit  assigned  for  tape 
containing  the  self-consistent-field  data 

NEDGE 

I 

=  -1:  terminates  Program 
=  0:  normal  run 

>  0:  sets  photon  energy  to  NEDGE  th 
binding  energy 

IA 

I 

=  0:  normal  case 

>  0:  calculation  commences  after  IA  shells 

IB 

I 

=  0:  calculation  to  Include  outermost 
shell 

>  0:  calculation  cuts  off  after  IB  shells 

IZ 

I 

Atomic  number 

L00P 

I 

=  0:  normal  case 

=  1:  photon  angular  momentum  reduction 

QV 

R 

Photon  energy  in  keV 

SAVE 

R 

-  3 

Total  cross  section  accumulation  from 
previous  run;  =  0  ordinarily. 

Input  Testing 


The  program  sifts  the  input  data  to  insure  that  certain  criteria  are 
not  violated.  If  any  difficulty  is  observed,  a  violation  signal  is  printed 
and  the  run  terminated.  Listed  below  in  abbreviated  form  are  the  criteria 
that  must  be  satisfied  and  that  the  program  tests  for: 

MIN  (JM,  KMAX,  IM,  IA,  IB)  2  0 
2  <  IZ  <  102 
QV  >  0 

L00P  =  0  or  1 
NEDGE  *  JX 
IA  <  JX 
IB  £  JX 
NK,  NKP  £  200 

In  addition  to  the  above  the  program  examines  the  following 
variables  and  alters  them  if  necessary: 

JM  >  24,  program  sets  JM  =  24,  prints  this  fact  and  proceeds. 

KMAX  >  12,  program  sets  KMAX  «  12,  prints  this  fact  and  proceeds. 

IM  >  12,  program  sets  IM  =  12,  prints  this  fact  and  proceeds. 

This  is  to  insure  that  the  dimension  size  of  the  variable  of  the 
program  as  written  for  the  IBM  7094  is  not  exceeded. 

Program  Options 

The  program  has  two  option  procedures  controlled  by  L00P  and  by  the 
combination  of  values  of  IA,  IB,  and  SAVE. 

By  setting  L00P  =  1,  we  can  reduce  the  original  angular  momentum 
quantity  LM  (max  J  for  photon)  in  unit  increments  (until  a  minimum  of  one  is 
reached).  The  selection  rules  may  then  reduce  the  range  limits  KM  and  JM.  The 
summing  process  is  repeated  to  recalculate  the  differential  and  total  cross 
sections . 
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By  altering  the  input  values  IA  and  IB  we  can  segment  a  run.  The 
program  proceeds  through  IB  electron  shells  rather  than  t.he  totality  of  shells 
(JX).  The  total  cross  section  up  to  that  point  can  be  led  back  in  as  part  of 
the  input  (SAVE)  when  the  user  desires  to  continue  the  run.  With  IA  greater 
than  zero,  the  program  starts  after  the  first  IA  electron  shells.  This  process 
is  useful  when  long  computer  running  time  is  not  available. 

System  Information 

The  program  as  run  at  the  New  York  Regional  Computer  Center  utilizes 
the  following  tapes: 

Logical  tape  unit  5  Read 

Logical  tape  unit  6  Write 

Logical  tape  unit  1  Self-consistent -field  data 
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III.  SAMPLE  OOTRJT 

The  output  from  an  actual  run  is  given.  The  complete  output  contains 
the  differential  cross  section  for  each  subshell  in  succession.  Only  one 
subshell  is  reproduced  here,  for  economy. 
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INPUT  data  ton  photoelectric  calculation 

1NHUI  CAKU  KtAUb 


UM  KMAX  Lfc  NIAHE  NEUGE  I A 
11  °  4  1  0  0 


Id 

0 


IZ 

13 


LOOP  ov 

0  0.10000000E  C2  G. 


SAVE 


INPUT  UA1 A  FOR  PHOTOELECTRIC  CALCULATION 


13  =  NUCLEAR  CHARGE 
1U.UUU  =  PHOTON  ENERGY 

o  =  MAX  KAPPA  FOR  ELECTRON 
4  =  MAX  L  FOR  PHOTO.*. 

11  =  MAX  J  (LEGENDRE  COEP'F.) 


TAPE  POSITIONED  PROPERLY. 


AL 

13 

6 

0.1944J3/0L  OU 
2b7 
1*61 
1436 

u.d*E221b2L  04 


ELEMENT 

ATOMIC  NUfcbCH 

NUMBER  OF  ELECTRON  SHELLS 

SCREENING  FACTOR  OF  OUTERMuST  lOUNO  ^lECT« 
RADIAL  GRIC  UP  TO  X  =  1.0 
RALIml  GRID  UP  To  X  =  65. u 

Total  radial  grid 

OUTERMOST  RADIAL  »/ALUL 


SHELL 

IMEbRATION  CUT-OFF 

KUi-HER  OF  WAVE  FUHCTICr.  GRID  POINTS 

p.aa  kappa  for  this  shell 

NU.-.HER  OF  MATRIX  ELEMENTS  FOR  THIS  SHELL 


LtNGTH  UNITS  ARE  HBAR  /  MC  (  1  oOHP.  RADIUS  =  137  ) 


INTEbKAT ION 

STEP 

SIZE 

IS 

0  •  CO  /fc  12b 

UP  TO 

l.OCO 

ANItGKAt ION 

STEP 

SIZE 

IS 

U  •  lifboUUU 

UP  TO 

65.000 

INI tGKA 1  ION 

STEP 

SIZt 

IS 

i.OUUCCOU 

UP  TO 

?.1Z .  693 
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*C  cr 
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li  =  MAX  J  (Lc.bi.NPNt  COi_FF.) 
b  =  MAX  KAPHA  NON 
*♦  =  MAX  L  FOK  PHOT i?f- 


LtbtNPKc.  COtFF'ICltNTb  OF  CNObS  StCTIO*i 
J  U(  J) 

U  U.64bd274bE  U 2 

1  U.d4b4114yt£  cc. 

2  -U./bdU00b2E  U2 

d  -U.dd2l6vbbt  U c 

4  -0.82VUb/t4i-:  U1 

b  -U.lb*b2244t  U1 

b  -U.22d4bblbE  UU 

/  -U.24odbUb4t-Ul 

-u.2UobyybbE-0i: 

-C.il/dV2.ddt-lb 


ELEMENT 

ATOMIC  NUMBER 

SHELL 

At. 

13 

1S1/2 

BINDING  ENERGY 

PHOTON  ENERGY 

ELECTRON  KINETIC  E>'EP3Y 

1.54S93U9  KLV 

10.0000000  KEV 

P. **500691  KEV 

UNR0LAK1ZEU  CROSS  SECTION  (I7ARNS/STERADIAN) 


lrttTA 

COS  THETA 

CROSS  SLOT  I CM 

ANG.  DIST. 

0 

1.0000000 

0.1 0482765* -01 

0.0000778 

2 

0.9993906 

0.31837447E 

0*7 

0.0023632 

4 

0. 9975640 

0. 123871 J2E 

01 

0.0091947 

6 

0.9945219 

0.27615944E 

01 

0.0204966 

b 

0.9902bbl 

0.487065635 

U1 

0.0361540 

10 

u. 9843078 

0.754335o5E 

01 

0.0559929 

12 

u. 9781476 

0.107513U3E 

02 

0.0798049 

14 

0.9702957 

0.14460679E 

C2 

0.10733P9 

lb 

0.9612bl7 

0.18632747E 

02 

0.1363074 

18 

0.9510565 

0.23224424E 

02 

0.1723906 

20 

0.939592b 

0.28188694E 

02 

0.2092409 

22 

u. 9271839 

0.33476231E 

Gc 

0.2484882 

24 

o .9135455 

0.39034332E 

02 

0.2897446 

2b 

0. 6987940 

0.44809127E 

G2 

0.3326098 

2b 

0.8629476 

0.507457375 

02 

0.3766766 

30 

U. 8660254 

0.56789170E 

0? 

0.4215355 

32 

0. 3480481 

0.628345alf 

02 

0.466780a 

34 

u. 829037= 

0.6897P249E 

02 

0.5120127 

3b 

t. 8090170 

0.7501P242E 

02 

0.5568465 

3b 

U. 7880108 

0.80954696E 

02 

0.6009117 

40 

u. 7660444 

0.6b7404olE 

G2 

0.6438584 

42 

o. 7431448 

0.92331445E 

0." 

0.6853596 

44 

0.7193396 

0.97667210E 

02 

0.7251140 

4b 

0.6946584 

0. 10277077E 

03 

0.7628484 

4b 

U .6691306 

0.10754931E 

03 

0.7933186 

SO 

0.6427676 

0.11199406E 

03 

0.831311 1 

S2 

U.olb6bl5 

0.1160P047E 

03 

0.e61642b 

S4 

0.5877853 

0. 11978619E 

03 

0.8891655 

bo 

0.5591929 

C.123101U7E 

03 

0.9137564 

S8 

0.5299193 

0.126007U2F 

03 

0.9353266 

bO 

U.  5000000 

0.12849794E 

03 

0.9530164 

62 

d  .4694716 

0.13056947E 

03 

0.969193C 

64 

0.4383712 

0 • 13222084E 

G^ 

0.9814503 

6b 

d  .40o7366 

0 .13345*. 

03 

0.9906067 

bb 

0 . 3746066 

0 . 13427632F 

Os. 

0.9967063 

70 

U. 3420202 

0.134d°443F 

05 

0.9998116 

n 

0.3090170 

0.13471976F 

03 

1.0000000 

74 

0.2756374 

0 .134365501 

C3 

0.9973702 

7b 

0.2419219 

0.13364bO0F 

Oj 

0.9920340 

78 

0. 2079117 

0.13257975F 

Os. 

0.9841149 

60 

U  •  1 73b452 

0.13118296E 

o: 

0.9737468 

82 

0.1391731 

0.129475J6F 

03 

C. 9610716 

84 

o. 1045285 

0.12747o93E 

0.9462376 
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"■ft**  '  VdCJ^Ti^pp^rj.’rj!  -— ■ 


H 

0*0697565 

0.12520626E  03 

0.9293979 

M 

0.0346995 

0.12269040C  D3 

0.9107082 

VO 

0.0000000 

0.11994453E  03 

0.8903260 

9t 

•0.0348995 

O.I1699190E  03 

0.8684091 

W 

-0.0697564 

0.11365362E  03 

0.8451144 

96 

-0.1045264 

O.I1055058E  03 

0.8205965 

96 

-0.1391731 

0. 1071032 3E  03 

0.7950075 

100 

-0.1736462 

0.10353156E  03 

0.7884956 

102 

-0.2079117 

0.99654979E  02 

0.7412050 

109 

-0.2419219 

0.96092227E  02 

0.7132748 

106 

-0.2756373 

0.92261351E  02 

0.6848389 

ioa 

-0.3090170 

0.88379639E  02 

0.6560257 

110 

-0. 3420201 

0.B4463592E  02 

0.6269576 

112 

-0.3746066 

0.8052889 IE  02 

0.5977510 

114 

-0.4067366 

0.765904U3E  02 

0.5685164 

116 

-0.4383711 

0.72662151E  02 

0.5393577 

118 

-0.4694715 

0.68757383E  02 

0 .5103729 

120 

-0.5000000 

0.64888321E  02 

0.4816540 

122 

-U. 5299192 

0.61066&79E  02 

0.4532866 

124 

-0.5591929 

0.57303176E  02 

0.4253509 

12b 

-0.5877852 

0.536078u5E  02 

0.3979206 

128 

-3.6156615 

0.49939820E  02 

0.2710652 

13U 

-0.6427876 

0.46457748E  02 

0.3448472 

132 

-0.6691306 

0.43019425E  02 

0.3193252 

134 

-0.6946533 

0.39682028E  02 

0.2945523 

136 

-0.7193398 

0.36452100E  02 

0.2705772 

138 

-0.7431448 

0.33335594E  02 

0.2474429 

140 

-0.7660444 

0.30337900E  02 

0.2251926 

142 

-C. 7880107 

0.274630/1E  02 

0.^038592 

144 

-0.8090170 

0 • 24717871E  02 

n. 1334762 

146 

-0.  <1290376 

0.221036C3E  02 

0.1640724 

146 

-0 .8480481 

0.19625136E  02 

0.1456737 

150 

-0.8660254 

0.17284947E  02 

0.1282030 

162 

-0.8829476 

0.15085945E  02 

0.1119602 

164 

-0.8987940 

0.130305GRE  02 

0.0967230 

186 

-0.9135454 

0 . 1112071  OF  02 

0.082547," 

188 

-0.9271638 

0.935834U0E  01 

0 . U694652 

160 

-U. 9396926 

0.77449520E  01 

C. 0574893 

162 

-0.9810565 

0.62B18ot<6E  01 

0.0466291 

lb4 

-C. 9612617 

0.49702C64E  01 

0.0368929 

166 

-U. 9702957 

0.3B109176E  01 

0. 0282377 

168 

-0.9781476 

0.2B047778E  01 

0.0208197 

170 

-0.9848078 

0.19524343F  01 

0.0144926 

172 

-0.9902681 

0.12544079E  01 

0.0093112 

174 

-0.9945219 

0.71109624E  00 

0 • 00627S3 

176 

-0.9975640 

0.32280832E  00 

0.0023961 

178 

-0. 9993908 

0.89763443E-01 

0 . 0006662 

180 

-1.0000000 

0.12053592E-01 

0 . 00003°5 

INTEGRATED  CROSS  SECTION  =  0.1066C397E  04  mar 
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ELEMENT  ATOMIC  NUMBER 

AL  13 

PHOTON  ENERGY  =  10.0000000  KEV 

TOTAL  CROSS  SECTION  =  0.11461591E  04  BARNS 

i 

i 


IV.  THE  PROGRAM 


The  program  is  written  in  FORTRAN  IV.  The  versions  used  on  the 
IBM  7091*  and  the  CDC  6600  are  identical  except  for  control  cards.  The  data 
tapes  contain  the  same  information  in  binary  form,  but  are  not  compatible 
between  the  two  machines. 

All  the  Common  blocks  appear  in  the  main  routine  (PELEC).  The 
definition  of  the  Common  variables  is  given  first.  In  the  subroutines.  Common 
blocks  used  are  quoted.  The  unlabelled  Common  is  the  same  wherever  it  appears. 


Definition  of  Variables  in  COMON 


Unlabelled  Common: 

Length 

iOTg 

Name 

Dimensions 

Mode 

Meaning 

PI 

R 

TT 

HA.LFPI 

R 

n/2 

FtfERPI 

R 

4  TT 

BAD 

R 

TT/l30 

SQ2 

R 

2'1/2 

Q 

R 

Photon  energy  (in  me  units) 

ZA 

R 

Atomic  number 

ZAZA 

R 

ZA  *  ZA 

EFN 

R 

2 

Free  electron  energy  -1  (in  me  units) 

EGN 

R 

2 

Free  electron  energy  +1  (in  me  units) 

V 

R 

Potential;  screening  factor/radius 

CG 

30 

R 

*k  -  Y  if  radius  <1;  -k  if  radius  >  1 

GAM 

30 

R 

'Ik2  -  ZAZA 

/BESSEL/  Common: 

Length 

WOg 

Name 

Dimensions 

Mode 

Meaning 

FL 

15 

"1 

)  Numerical  factors  used  in  the 

►  construction  of  the  spherical 

PC 

15 

R 

!  Bessel  function. 

0F 

15  x  15 

R  * 

1 

Ml 

I 

Largest  order  of  Bessel  function  needed 

M2 

I 

Ml  +  1 

B 

15 

R 

Spherical  Bessel  function 

/DFUNC/  Cannon:  Length  1047g 


Raine 

Dimensions 

Mode 

Meaning 

F 

30 

R 

"Small"  component  of  free-electron 
wave  function 

G 

30 

R 

"Large"  component  of  free-electron 
wave function 

0F 

30 

R 

Derivative  of  F 

DG 

30 

R 

Derivative  of  G 

DFK 

200 

R 

Integrand  for  matrix  elements  K  (kx') 

X 

DFKP 

200 

E 

Integrand  for  matrix  elements  K  (h’h) 

X 

CF 

30 

R 

h  -  y  if  radius  <  1;  h  if  radius  >  1 

H 

R 

Integration  step  size 

/FAC /Common 

• 

• 

Length 

32kg 

Name 

Dimensions 

Mode 

Meaning 

FACT 

67 

R  “ 

Numerical  factors  used  in  the  calcula¬ 
tion  of  the  Clebsch-Gordan  coefficients 

KTFACT 

95 

R 

R00T 

50 

R  - 

/FID0/  Common: 

Length 

7k6g 

Name 

Dimensions 

Moc, 

Meaning 

FI 

30  x  15 

R 

0(k,\) 

D 

30 

R 

Legendre  coefficient  of  cross  section 

JMP 

I 

Max.  order  of  Legendre  coefficient  +  1 

NAME 

ALFA 

Element 

SHELL 

ALFA 

Electron  shell 

QV 

R 

Photon  energy  (in  keV  units) 

EB 

R 

Binding  energy  of  shell  (in  keV  units) 

IZ 

I 

Atomic  number 
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/HUT/  Cannon: 


Length  24g 
Mode  Mean: 


Name 


Dimensions 


ing_ 


RK1 

4 

R  "1  Numerical  coefficients  used  for 

■"  the  Runge-Kutta  integration 

RK2 

4 

R 

(Gill  Form) 

RK3 

4 

R 

RK4 

4 

R 

K4 

4 

I  - 

/LIMIT/  Cannon: 

Length 

1J6 

Name 

Dimensions 

Mode 

Meaning 

JM 

I 

Max  order  of  Legendre  coefficient 

IM 

I 

Max  j l  for  photon 

KM 

I 

Max  k  for  free  electron 

K2M 

I 

2  *  KM;  number  of  free  electron  states 

IEND 

I 

Flag  for  zero  electron  kinetic  energy  state 

NEW 

I 

Flag  to  save  repetition  of  calculation  when 
radius  is  not  advanced 

NK 

I 

Number  of  matrix  elements  K  (hh’) 

At 

NKP 

I 

Number  of  matrix  elements  K  (k'k) 

JKB 

I 

Twice  Jh, 

IMKB 

I 

NTAB 

I 

Radial  index  in  bound  state  tabulation 

/MAT/  Common: 

Length 

15358 

Name  Dimensions 

Mode 

Meaning 

SF 

30 

R 

Integration  storage  variable  for  "small" 
component  of  free  electron  wavefunction 

SG 

30 

R 

Ingegration  storage  variable  for  "large" 
component  of  free  electron  wavefunction 

FK 

200 

R 

Matrix  elements  K.(kh') 

FKP 

200 

R 

Matrix  elements  K.(k'k) 
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Rase 

Dimensions 

Mode 

Meaning 

SFK 

200 

R 

Integration  storage  variable  for  matrix 

elements  K  (hh') 

fL 

SFi-c* 

20C 

R 

Integration  storage  variable  for  matrix 
elements  K.(h'h) 

RCUT 

R 

Cut-off  radius  for  integration 

/QUART/  Cannon: 

Length  322g 

Name 

Dimensions 

Mode 

Meaning 

LK 

50 

I 

IKK 

50 

I 

'-H 

JK 

50 

I 

Twice 

FKAP 

50 

R 

H 

SN 

50 

R 

Sign  of  k 

SI 

50 

R 

Sine  of  phase  shift 

CR 

50 

R 

Cosine  of  phase  shift 

/ONWARD/  Common: 

Length 

U8 

Name 

Dimensions 

Mode 

Meaning 

RX 

R 

Radial  variable  beyond  radial  cutoff  on 
uniform  grid  tabulation 

sex 

R 

Interpolated  value  of  screening  factors 

GBX 

R 

Interpolated  value  of  "large"  component 
of  wave function 

FBX 

R 

interpolated  value  of  "small"  component 
of  wave function 

/TAPES/  Common: 

Length 

Name 

Dimension 

Mode 

Meaning 

X 

1500 

R 

Radial  value 

SCF 

1500 

R 

Screening  factor 

F3 

1500 

R 

"Small"  component  of  bound  state 

wave function 
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Bane 

Dimensions 

Mode 

Meaning 

GB 

1500 

R 

"Large"  component  of  bound  state 
wave function 

GAMB 

R 

k2  -  ZAZA,  for  bound  state 

SCREEN 

R 

Normalization  factor  for  least  bound 
electron 

/TRANS/  Ccaanon: 

Length 

1700g 

Name 

Dimensions 

Mode 

Meaning 

HF 

30  x  15 

R 

H  (h  ,(! 1 )  for  n*  >0 

HFM 

30  x  15 

R 

H  (h j|i * )  for  (i*  <  0 

JNG 

30 

I 

(-p‘)  max  +  l/2 

JPS 

30 

I 

(u*)  max  +  1/2  provided  p'  >  0  permitted 
-1  otherwise 

/VECT/  Common: 

Length 

Mh 

Name 

Dimensions 

Mode 

Meaning 

KF 

200 

I 

Index  for  a  values  for  K^(h'h) 

KG 

200 

I 

Index  for  h  values  for  K  (h  h') 

At 

LBES 

200 

I 

Photon  angular  momentum  +  1  for  K.(h’  h) 

At 

LBS 

200 

I 

Photon  singular  momentum  +  1  for  K .  (h  h  1 ) 

At 

LKB 

I 

l  , 
h’ 
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ROUTINE  FEXEC 


Purpose:  This  is  the  main  program.  It  serves  as  the  control  section, 

reading  the  input  cards  and  the  tape  containing  the  self- 
consistent-field  data  and  setting  up  the  calculations.  Nearly 
all  the  actual  computation  is  done  in  subroutines. 

Method:  Numerical  constants  are  computed  and  stored,  the  input  read  and 

checked  (and  diagnostics  printed,  if  needed),  then  printed  with 
interpretation  the  tape  is  scanned  for  the  element  required  and 
the  data  read  from  it  and  interpreted.  For  each  subshell, 

SINDEX  is  called  to  set  up  the  matrix  elements  and  their  quantum 
numbers,  RADINT  to  perform  the  radial  integrations.  LEGEND  is 
called  to  perform  the  angular  momentum  sums,  ANGLE  to  produce 
the  differential  cross  section  (headings  are  supplied  for  the 
tables  printed  in  the  subroutines);  if  the  matrix  reduction 
option  is  called  for,  the  cutoff  values  of  the  quantum  numbers 
are  reduced  and  the  program  loops  back  through  the  angular  stuns. 
The  total  cross  section  is  accumulated  in  PELEC  and  printed 
out.  The  program  recycles  to  read  the  next  input  card  and 
perform  the  next  case,  until  a  flag  on  the  input  card  signals 
termination. 

Subroutines  called:  SINDEX,  RADINT,  LEGEND,  ANGLE,  HUM 

Variables  in  unlabelled  Common:  PI,  HALFPI,  F0URPI,  RAD,  SQg,  Q,  ZA,  ZAZA, 

EFN,  EGN,  V,  CC-,  GAM 

Labelled  Common:  BESSEL,  DFUNC,  FAC,  FUtf,  KUT,  LIMIT,  MAT,  QUANT,  INWARD, 
TAPES,  TRANS.  VECT 

3  * 


Local  Variables 


Name 

Dimension 

Mode 

Meaning 

KMAX 

I 

Input  max  *t  for  free  electron 

NTAPE 

I 

Hie  logical  tape  unit  assigned  for  tape 
containing  the  self-consistent-field 
data 

HEDGE 

I 

=  -1  terminates  program 
=  0  normal  run 

=  integer  sets  photon  energy  to  NEDGEth 

binding  energy 

IA 

I 

-  0  normal  case 

>  0  calculation  conmences  after  IA  shells 

IB 

I 

=  0  calculation  includes  outermost  shell 
^  0  cuts  off  calculation  after  IB  shells 

Lfi»P 

I 

=  0  normal  case 

=  1  for  photon  angular  momentum  reduction 

QV 

R 

Input  photon  energy  in  keV 

SAVE 

R 

Input  total  cross  section;  accumulation 
from  previous  run;  =  0  otherwise 

IZ 

I 

Input  atomic  number 

lgrid 

I 

Lumber  of  grid  points  in  table 

ISDCFI 

I 

Number  of  grid  points  for  radial  value 
of  (l/2  Bohr  unit) 

ItfUE 

I 

Number  of  grid  points  for  radial  value  of 
one 

RSIXFI 

R 

Radial  value  of  6 %  read  from  tape 

ROHE 

R 

Radial  value  of  one,  read  from  tape 

ZTRY 

R 

Atomic  number  read  from  tape 

JX 

I 

Number  of  electron  shells 

XN 

36 

ALFA 

Shell  identification 

XL 

36 

R 

iH> 

XJ 

36 

R 

V 

XZ 

36 

R 

Shell  occupancy 

Heme 


Dimension 


Mode 


^  ERG  36  R 

KJI  I 

SECT  R 

NT0T  I 

SEDGE  R 

IMS  I 

JMS  I 

JPS  I 

CCM  R 

ALFA  R 

REL  R 

KB  I 

C  FKB  R 


Meaning _ 

Binding  energy 

Number  of  wavefunction  grid  points 

Total  cross  section 

Number  of  matrix  elements  for  shell 

Cross  section  jump  at  edge 

Initial  max  t  (in  loop  reduction) 

Initial  max  order  of  Legendre  coefficient 
(in  loop  reduction) 

Initial  max  order  of  Legendre  coefficient 
+  1  (in  loop  reduction) 

2 

Conversion  factor  (keV/mc  ) 

Fine  structure  constant:  1/137.0367 

12 

Classical  electron  radius  (in  cm)  x  10 

V 

V 
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FELBC  > 
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r PELEC 
Page  4  of  6 


Energy  falls  below  edge; 
bypass  computation 


Is  photon  energy  =  binding 
energy 


Set  up  quantum  numbers,  energies, 
cutoff  values  for  shell 


Is  computed  max  h  for  free  electron  = 
input  value  of  h 


Write  out  computed  value  which  will  be 
the  one  used 


Set  NEDGE  *  loop  index 


\  » 

Calculate  no.  of  free 

/  * 

electron  states 

fcall  subroutine  SHJDE» 

V  _ 

itrix  element  dimensions  exceeded 
|  No 

Commence  integration \ 
Call  subroutine  RADINT 

1 

*\^Cal2subrout  ine  LEGEND  \ 


Multiply  Legendre  coef  by 
numerical  factor 


Call  subroutine  ANGIE 


Write  out  this 
information 


Loop  on  max  order 
'of  leg.  coef  +  1 
to(o6) 


c 


tern  reduction  flag 
<  2 


i 


No 


Is  max  i  -  1 


pi— © 


I 


No 


Reduce  max  t  by  one 


I 


(Is  max  h  for  free  electron  £ 
(h'  +  max  t) 


No 


Reset  max  k  no.  of  free  electron 
states ,  and  max  Legendre  coefficient 


< 


i 


Start  summing  process 
Call  subroutine  HUM 


© 


> 
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w 


Set  flag  -  1 

reset  previously  stored  max  l , 
max  Leg.  coef 


o  r 


51BMC  PELEC 

►OH  CUC  6600  OPERATION.  ADO  CARD  REA0XN6 

PR06KAK  PELEC ( INPUT  .  OUTPUT  #  TAPE1 *  T APE 5= INPUT . TAPE6=0UTPUT ) 

COMMON  PI  .FIALFPI.FOUHPl»KAD.SQ2»Q.ZA.ZA2A.EFN»E6N»V.CG(30) »GAM(30)PEL00010 
COMMON/BEbbEL/FL( lb) »PC(15) .0F(15»15) »M1.M2.B(15)  PEL00020 

COMMON/DFUNC/F(30)  »G(30)  .OF (30)  *0600)  »OFK( 200)  »0FKP(200 ) .CF( 30)  'HPEL00030 


COMMON  /FAC/FACT (67) 'RTF AC (9b)' ROOT (50) 

COMMON/ F IUO/F 1 ( 30 » 1 b ). U ( 30 ) t JMP . NAME » SHELL » QV » EB » IZ 
COMMON  /KUT/RKl (4) .RK2(4) .RK3(4) »RK4(4) »K4(4) 

COMMON  /LIMIT/UM.LM.KM.K2M. IEND.NEW.NK.NKP. JKB’LMKB.NTAB 


PEL00040 

PEL00050 

PEL00060 

PEL00070 


COMMON/MAT/SF (30) .S6(30) »FK(200) »FKP(200) »SFK(200) »SFKP(200) .RCUT  PEL00080 


COMMON/QUANT /LK (30) »LMK(30) » JK(30) »FKAP(30) »SN(30) .51 (30) .CRT 30) 
COMMON  /ONMARD/KX.bCX.GbX.FBX 

COMMON  / TAPtS/X ( lbOO ) » SCF ( 1500 ) »FB( 1500 ) .00(1500 ) . 6AM0 . SCREEN 
COMMON/ TR ANb/hF (30.15). HFM (30.15). JNG ( 30 ) . JPS ( 30 ) 

COMMON/ VEC I /KF (200) »Kti(200) . LBES(200) » LBS (200) .LKB 
UIMtNblON  XN ( 36 ) . XL ( 36 ) . X J ( 36 ) *ER6(36) .XZ(36) 

FL(1)  =  1.0 

PC(1)  -  1.0 

UO  5  L=2»l5 

FL(L)  =  2*L-1 

PC(L)  =  PC(L-1)/FL(L) 

UO  5  J=1*15 

FLU  =  J*(2*(L+J)-1) 

5  OF(L.J)  =  l.O/FLJ 
KOU 1(1)  =  1*0 
UO  1U  1=2.50 
FAT  =  1 

1U  KOOI (1)  =  50RT  (FAT) 

FACI (1)  =  1.0 
F  AC  I ( 3 )  =  1.0 
KTFAC(l)  =  l.U 
KTF AC (3)  =  l.U 
FAT  =  l.U 
UO  15  1=2.33 
FI  =  1 

FAT  =  F A  I *F I 
FACI (2*1+1)  =  FAT 

15  HTFAC(2*I+1)  =  bQRT(FAT) 

FAT  =  l.U 
UO  2U  1=34.47 
Fi  =  I 

FAT  =  FATAFI 

2U  KTFAC(2*I+1)  =  bQRT (FAT ) *RTFAC (67) 

502  =  1.0  /  MOOT (2) 

KK1  (1)  =  U.5 
KK1  (2)  =  1.0-502 
KK1  (3)  =  1.0+502 
RK1  (4)  =  1.0/6. 0 


=  1.U-SU2 
=  1.0+502 
=  U.5 
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K4(l)  = 


PEL00090 
PEL00100 
PEL00110 
PEL00120 
PEL00130 
PEL00140 
PEL00150 
PEL 00 160 
PEL00170 
PEL00180 
PEL00190 
PEL00200 
PEL00210 
PEL00220 
PEL00230 
PEL00240 
PEL00250 
PEL00260 
PEL00270 
PEL00280 
PEL00290 
PEL00300 
PEL00310 
PEL00320 
PEL00330 
PEL00340 
PEL00350 
PEL00360 
PEL00370 
PEL00380 
PEL00390 
PEL00400 
PEL00410 
PEL00420 
PEL00430 
PEL00440 
PEL00450 
PEL00460 
PEL00470 
PEL00480 
PEL00490 
PEL00500 
PEL00510 
PEL00520 
PEL00530 
PEL00540 
PEL00550 
PEL00560 
PEL00570 
PEL00580 
PEL00590 


K4C2)  =  0 

K4I3)  =  1 

*414}  =  U 

CCM  =  SU.UUbZ 

ALFA  =  1.U/137.U367 

KtL  =  U. 281 777 

PI  S  3.141S926S 

FOUHPI  =  4 . U*PI 

HALF P I  =  U.b  •  PI 

KAU  =  PI  /  18U.U 

FIB  =  U.2bU*PI*REL*HtL/ALFA 

22  Hk.AU  (S#2S)  jm # kmax # lm#nt ape » nedge#  IA'  ib»  iz » loop #ov# save 

Zb  FOHMAI  (91b'2tlb.8) 

C  OHU1NAHILY t  IA  =  0*  lb  =  0#  SAVE  =  0. 

C  lb  .Nt.  U  CUTS  OFF  CALCULATION  AFTER  IB  SHELLS 

C  1A  .01.  U  RtSTAHTS  AFTER  IA  SHELLS 

C  SAVE  ib  IHLN  CROSS  SECTION  OF  FIRST  I A  SHELLS  READ  BACK  IN 

C  uv  ENIEKEU  IN  K.E.V.I  SAVE  ENTERED  IN  BARNS 
C  LOOP  =  0  ONUINANILY »  LOOP  =  1  FOR  LM  REDUCTION 
C  NtUOE  NEGATIVE  TERMINATES  PROGRAM 
IF  (NtUGt.LT.U)  CALL  EXIT 

WK1IE  (b»2/) 

2/  FORMAT  (1H1///2X.4UHINPUT  DATA  FOR  PHOTOELECTRIC  CALCULATION 

LM  NTAPE  NEDGE  IA  IB 


PEL00600 

PEL00610 

PEL00620 

PEL00630 

PEL00640 

PEL00650 

PEL00660 

PEL00670 

PEL00680 

PEL00690 

PEL00700 

PEL00710 

PEL00720 

PEL00730 

PEL00740 

PEL00750 

PEL007b0 

PEL00770 

PEL00780 

PEL00790 

PEL00800 

PEL00810 

//2X#PEL00820 


HbHlNPUT  CARD  KEAUS//2X * 38HJM  KMAX 
2bH  LOOP' /X'2HUV» 1UX.4HSAVE  //) 

WRITE  (b»2b)  JM  t KMAX  >  LM » NT  APE » NEDGE  'lAtlB'IZ*  LOOP  » QV • SAVE 
UAIA  CHECKING 

IF  (MlNU(JM'KMAXfLM.IA'lB).GE.U)  GO  TO  29 

2b  WRITE  (b» 1U2) 

1U2  F OHM AT  (bX» lbHlNCOHRECT  INPUT  //I 
GO  10  22 

29  IF  ((1Z.LI.2)  .OH.  (ll.GT.102)  )  GO  TO  36 

(JV.Lt.U.U)  GO  TO  28 

( (LOOP. NE.U).ANU. (LOOP. NE.l))  GO  TO  28 


1U9 


GO  10  111 


IF 
IF 
IF 

IF  (JM.Lt.24) 

JM— 24 

NKlIt  (b» 109)  JM 

FORMAT  (/bXbbMlN  ORDER  NOT  TO  EXCEED  DIMENSION  JM  HAS  BEEN 
1  10  13/) 

111  IF  (KMAX.Lt.12)  GO  TO  112 

KMAX— 12 

WHIlt  (b» 113)  KMAX 

113  FORMAT  (/SXS8H1N  ORDER  NOT  TO  EXCEED  DIMENSION  KMAX  HAS 
1EU  10  13/) 

112  IF  (LM.LL.12)  GO  TO  116 

LM— 12 

WHIlt  (b» 114)  LM 

114  FORMA  I  (/bXbbHIN  ORDER  NOT  TO  EXCEED  DIMENSION  LM  HAS 
1  TO  13/) 

lib  2  =  12 

2A  -  2*ALF A 
ZAZA  =  ZA*ZA 

WRITE  IbfbU) IZ'QV'KMAX' LM' JM 

bU  FORMAT  (1H1///SX40H1NPUT  DATA  FOR  PHOTOELECTRIC  CALCULATION  // 

1 14X f  14 » 1 7h  =  NUCLEAR  CHARGE  /9X»F9.3»16H  =  PHOTON  ENERGY  // 
21SX13»2SH  =  MAX  KAPPA  FOR  ELECTRON  /15XI3'19H  =  MAX  L  FOR  PHOTON 
31SX13'26H  =  MAX  J  (LtGENURL  COEFF.)  //) 

REW1NU  NTAPE 

READ  (NTAPE)  NOR  ID*  1SIXF  1 » IONL  'RSI  XF I  'RONE  r  ( X  ( I ) ' 1  =  1 ' hlGRID) 

3b  READ  (NTAPt)  /TRY » JX » NAME t SCRtEN' ( XN( I ) ' XL ( I ) ' X J( I ) » XZ ( I ) ' ERG( I ) ' 
11=1' JXJ 

12 1  =  ZTRT  ♦  U.U1 


IZ»PEL00830 
PEL00840 
PEL008S0 
PEL00860 
PEL00870 
PEL00880 
PEL00890 
PEL00900 
PEL00910 
PEL00920 
PEL00930 
PEL00940 
PEL009b0 
PEL00960 
REDUCEDPEL00970 
PEL00980 
PEL00990 
PEL01000 
PEL01010 
BEEN  REDUCPEL01020 
PEL01030 
PEL0104U 
PEL01050 
PELOlObO 
BEEN  REDUCEDPEL01070 
PELOlObO 
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PEL01090 

PEL01100 

PEL01110 

PEL01120 

PEL01130 

PEL01140 

/PEL01150 

PEL01160 

PEL01170 

PEL01180 

PEL01190 

PEL01200 

PEL01210 


Gb(l)»l=l.KJI  *»  (  FBI I  * » 
RECORDS  TO  NEXT  ELEMENT 
GO  TO  35 


I=1»KJI  ) 


IF  l  U  .tu.  1ZT  *  60  TO  40 

KLAU  INIAPt*  ISCF ( X ) . 1=1 .N6KIU* 

SKIPS  LOGICAL  HECOHU  CONTAINING  POTENTIAL 
00  99  1SK1P=1.JX 

99  Kfc.AU  INTAPt*  RCUT.KJI.l 
SKIPS  UNWANTED  LOGICAL 
IF  (  UT  .Lb.  102  ) 

3b  NKlIt  lb»3U*  U 

30  FOKMAl  l//3UX19HNUCLtAK  CHARGE  Z  =  I3.1X18HN0T  FOUNO  ON  TAPE.* 
GO  10  22 

40  NKlIt  lb»45* 

4b  FORMA I  l/////2UX.25HTAPt  POSIIIONLO  PROPERLY.////) 

IF  (NtUGt. GT.UX*  GO  TO  26 

IF  (  UA.Gt.JX*  .OH.IIB.GT.JX*  *  GO  TO  28 

KLAU  (NTAPb)  ISCFll* »  l=l.NGKiU) 

IF  (  NtUGt  .bO.  0  J  GO  TO  93 
0  =  tHG( NtUGt* 

UV  =  U  *  COM 
GO  10  94 

93  U  =  UV  /  CLM 

94  WRilt  lb. 1/3)  NAME# 1/ »JX . SCREEN# IONE.RONE. ISIXFI .RSIXFI »NGRIO» 
IX (NGKlU* 

1/3  FORMA  I  1//25XA6»1X9H=  ELtMENT/28XI3. 1X1SH=  ATOMIC  NUMBER/ 
129X12. 1X2 7h-  NUMBtR  OF  ELECTRON  SHELLS/ 

216XtlS«8. 1X46H—  SCREENING  FACTOR  OF  OUTERMOST  BOUND  ELECTRON/ 
328X13. 1X24H=  RADIAL  GRID  UP  TO  X  =  F4.1/ 

RADIAL  GRID  UP  TO 
TOTAL  RADIAL  GRID 


X  =  F  b. 1/ 

/16XE15.8.1X24H=  OUTERMOST  RADIAL 


IB  =  JX 


427X14»1X24H= 

52  /  A 14. 1X19H— 

SALUt  ✓/) 
btCI  =  SAVt 
IF  l  lb  .tU. 

DO  bU  1  =1. 

REAU  INTAPt) 

IF  U.Lt.lA) 

IF  lU.Gt.tRGll) ) 

WRilt  lb»95)  XN ( I ) 

9b  FORMAT  l/bX.lBHENtKGY  FALLS  BELOW . 1 X A6 . 1 X4HEUGE// ) 
GO  10  bU 

9U  IF  l  Q  .tO.  ERG l 1 )  )  NtUGt  =  1 
*R11E  lb. 80*  XNIII .HCUT.KJI 


0  ) 

IB 

RCUT.KJl.  (GB(L) »L=1.KJI) . (FBIL* »L=1»KJI* 
GO  TO  bU 

GO  TO  90 


SHELL/16XE1S.8»1X21H=  INTEGRATION  CUT-OFF/ 
WAVE  FUNCTION  GRID  JOINTS  /> 


80  FORMAT llhl/2bXAb»lx7H= 

12 7X14. 1X3 /H=  NUMBER  OF 
FST  =  Flb*X2ll*/Q 
LKb  =  XLll)  ♦  0.01 
JKb  =  2 • 0*X J ( 1 )  ♦  0.01 
LMKb  =  JKB-LKB 
FKd  =  XJID+O.b 
GAMb  =  SOR I l FKB**2~2A*2A ) 

KB  =  FKb  ♦  0.01 
KM  =  MAXO  l  Kb+1»  KMAX  * 

KM  =  MINU  l  KM.  KB+LM  ) 

IF  IKM.EO.KMAX)  GO  10 

WRilt  l b » 18  7 )  KM 

187  FORMA  1  128X13. 27H  =  MAX  KAPPA  FOR  THIS  SHELL  ) 
18b  K2M  =  2  *  KM 

JMP  =  MINU  l  JM+1.K2M  ) 

EFN=U-ERGll ) 
td=tRG( 1 ) *LLM 
SHtLL-XNl 1 ) 

CALL  SINUtX 

IF  l  INK. Lt. 200*  .AND.  UmP.Lt.200*  ) 

WK 1  I b  Ib.iOU)  NK.NKP  .  _ 


18b 


GO  TO  83 


PEL01220 
PEL01230 
PEL01240 
PEL012S0 
PEL01260 
PELG1270 
PEL01280 
PEL01290 
PEL01300 
PEL01310 
PEL01320 
PEL01330 
PEL 01 340 
PEL01350 
PEL01360 
PEL01370 
PELC13oO 
PEL01390 
PEL01400 
PEL01410 
PEL01420 
PELOinjO 
PEL01440 
PEL01450 
PEL014b0 
PEL01470 
PEL01480 
VPEL01490 
PEL 01 So 0 
PEL01510 
PEL01520 
PEL01530 
PEL01540 
PEL01550 
PEL01560 
PELC1570 
P£L015o0 
PEL01590 
PELOlbOO 
PEL01610 
PEL01620 
PEL01630 
PEL01640 
PEL016b0 
PEL01660 
PEL01670 
PEL01680 
PEL01690 
PEL01700 
PEL01710 
PE'  "  ^0 

PELL  0 
PELOi <  40 
PEL017b0 
PEL01760 
PELOI 770 
PEL01780 
PEL01790 
PEL01800 
PEL01810 
PEL01820 
PEL01830 


. v  a*  i  y 


1UU  FORMAT  tSX39HMATRIX  ELEMENT  DIMENSIONS  EXCEEDED  *  1  4XSHNK  =  14*  PEL01840 

I6HNKP  =  I4//30X12HCASE  DROPPED  //)  PEL018S0 

w  10  «  PEL01860 

83  NTOl^NKFNKP  PEL01870 

•KITE  (6*1101)  NTOT  PEL01880 

1101  FORMAT  (27X»I4*43H  =  NUMBER  OF  MATRIX  ELEMENTS  FOR  THIS  SHELL  //>  PEL01890 
CALL  RAUINT  PEL01900 

88  CALL  LEGEND  PEL01910 

DO  86  J=1*JMP  PEL01920 

86  D( J)sFST*D( J)  PEL01930 

■RUE  (6*2  )  JM*  KM*  LM  PEL019**0 

i  FORMAT  UH1///1SX*I3*26H  =  MAX  J  (LEGENORE  COEFF.)/  15X*13*  PEL01950 

1  2SH  =  MAX  KAPPA  FOR  Li.ECTR0N/lSXl3*  19H  =  MAX  L  FOR  PHOTON/)  PEL01960 

CALL  ANGLE  PEL01970 

IF  (L00P.LT.2)  GO  TO  88  PEL01960 

•Kilt  (6*62)  LM  PEL01990 

62  FORMAT  (////S0X23H  LOOP  REDUCTION  TO  LM  =  12)  PEL02000 

6/  IF  (LM.EO.l)  GO  TO  63  PEL02010 

LM  =  LM-1  PEL02020 

IF  IKM.LE.lKb+LM) )  GO  TO  ol  PEL02030 

Ml  =  Kb+LM  PEL02040 

K2M  =  2*KM  PEL02050 

JMF  =  MINI)  (  JMP*  K2M  )  PEL02060 

DM  =  JMP  -  1  PEL0207C 

61  CALL  HUM  PEL02080 

GO  IJ  8b  PEL02C90 

88  IF  U.EG.NEDGE)  SEDGE  =  F0URPI*D(1)  PEL0210C 

SEC  I  =  SECr+FOURPI*D(l)  PEL02110 

IF  ILOOP.EU.U)  GO  TO  60  PEL02120 

LOOP  =  2  PEL02130 

CMS  =  LM  PEL02140 

JMS  =  JM  PEL 02 ISO 

JPS  =  JMP  PEL02160 

GO  10  87  PEL02170 

63  LOOP  =  1  PEL02160 

LM  =.  LMS  PEL02190 

JM  =  JMS  PEL02200 

JMP  =  JPS  PEL02210 

6U  CONIINUt  PEL02220 

«R1 lt(6»20U)  NAME*  12*  GV*  SECT  PEL02230 

2UU  FORMAmMl////30X*7HtLEM£NT»27X»13HAT0MIC  NUMBER//32X*  A6»31X»  I3///PEL02240 

1  30X.17HPHOTON  ENERGY  =  »F12.7»bH  KEV  ////  40X*  PEL02250 


2  22HI01AL  CROSS  SECTION  =  #E15.8*7H  BARNS  ////  ) 

IF  (NtUbt.tU.U)  GO  TO  22 

IF  A I  EDGE*  CROSS  SECTION  BELOW  EDGE 

SEC  I  =  SEC 1 -SEDGE 

WRICE  (6»/S)  SECT 

TS  FORMAI  (/4UX27HCR0SS  SECTION  uELOW  EDGE  =  E1S.8»6H  BARNS/) 
GO  10  22 
END 


PEL02260 

PEL02270 

PEL02280 

PEL02290 

PEL02300 

PEL02310 

PEL02320 

PEL02330 


-  33  - 
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SUBROUTINE  ANGLE 


Purpose:  This  routine  computes  the  angular  distribution  from  the 

Legendre  coefficients  of  the  cross  section  and  a  calculation 
of  the  Legendre  polynomials,  and  outputs  the  differential  cross- 
section  at  two  degree  intervals. 

Method:  Legendre  polynomials  are  computed  at  two  degree  intervals ,  then 

multiplied  by  the  Legendre  coefficients  of  the  cross  section  and 
finally  sunned. 

Subroutine  called:  None 

Subroutine  called  by:  PEIEC 

Variables  in  unlabelled  Common:  PI,  HALFPI,  PgURFl ,  RAD,  3G2,  ZA,  ZAZA, 

EFN,  V,  CG,  GAM 

Labelled  Common:  FID0 

Local  Variables: 


Name 

Dimension 

Mode 

Meaning 

JL 

I 

Order  of  Legendre  coefficient 

EL 

R 

Electron  kinetic  energy  in  keV 

SI-AX 

R 

Maximum  differential  cross  section 

KMUQ 

I 

Inde.xing  variable 

I THETA 

I 

Angle  (in  degrees)  between  photon  and 
electron 

IT 

91 

I 

Stores  angle  values 

EM uq 

R 

Cosine  of  angle 

SIP 

R 

Unpolarized  cross  section  (in  barns/ 
steradian) 

PPL 

R 

Legendre  polynomial 

3P 

91 

R 

Stored  values  of  unpolarized  cross  secti 

AD 

91 

R 

Angular  distribution 

XGRC 

R 

Integrated  cross  section  (in  barns) 

QMU 

91 

R 

Stored  values  of  cosine  of  angle 

35  - 

Write  Legendre  coefficients 
of  cross  section 


SltfMC  ANGL 

SUBROUTINE  ANGLE  ANGQ0010 

COMMON  P I . HALFP 1 . FOUKH I *  KAO  *  S02  *0  #  Z A «  ZAZA  *  EFN  #  EGN  » V  #C6 ( 30 ) • GAM ( 30  > ANG0002 0 


COMMON/F IOO/F 1( 30 » 1 5 ) . 0 ( 30 ) » JMP » NAME . SHELL * OV » EB » IZ 
DIMENSION  IT (91) *  GMU191) .  SP(91>*  A0(91) 

MRITt  (6*2  ) 

2  FORMAT (///1UX.38HLEGENDHE  COEFFICIENTS  OF  CROSS  SECTION// 

1  bX* 1HJ* 10X*4HD( J)  //  ) 

DO  4  J  r  1»  JMP 
JL  =  J  “  1 

WRIIE  lb>)  )  JL»  U(J) 

3  forma i  (  i7.sx*  eis.8  ) 

4  CONTINUE 

EL  =  GV  -  tb 

WRITE  (b» S)  NAME*  12*  SHELL*  LB*  GV*  EL 

5  FORMAT  (1H1//// »33X* 7HELEMENT » 7X* 13HAT0MIC  NUMBER*  13X*5HSHELL// 

1  3SX » Ao» 11X » 13* 17X* A6/ //lbX* 14H8IN0ING  ENERGY. 17X* 

2  13HPHOTON  ENERGY . 11X . 23HELECTRON  KINETIC  ENERGY// 

3  3( 14X *F 12* / »4H  KEV)  ///  ) 

IF  (  JMP  .tO.  1  )  GO  TO  24 

WRilt  (b*b> 

a  FORMAI (2SX.43HUNP0LARIZE0  CROSS  SECTION  ( BARNS/STER AD I AN )  // 

1  lbX  *  SMl HET  A* lbX  *  9HC0S  THETA *18X» 13HCR0SS  SECTION*  lbX* 

2  1UHANG.  01S1.  //  ) 

SMAX  =  O.U 

UO  IS  KMUG  =  1*  91 
1 THETA  =  2  *  (  KMUG  -  1  ) 

IT (KMUG)  =  1TMETA 
IHE=RAU*FLOAT  (I THETA) 

FMUO=COS  (I HE) 

GMU(KMUG)  =  FMUG 

PM1=1.0 

PN=FMUG 

SIP  =  Ull)  ♦  FMUG  *  U(2) 

IF  (JMP. EG. 2)  GO  TO  12 

UO  1U  J-3* JMP 

FUSJ-2 

PPL  =  (  PN*FMUG*(2.U*FN+1.U)  -  PMI*FN  )  /  (  FN+1.0  ) 

PM1-PN 

PN=PPL 

10  SIP  =  SIP  +  PPL  *  U(J) 

12  SP (KMUG )  =  SIP 

IF  (  SP(KMUG)  .GT.  SMAX  )  SMAX  =  SP(KMUQ) 

IS  CONIINUE 

UO  20  KMUQ  =  1*  91 
AD (KMUG)  =  SP(KMUQ)  /  SMAX 

WRITE  (  b.ltt)  IT (KMUG) *  QMU(KMUQ) »  SP(KMUQ).  AD (KMUQ) 


>  AD (KMUQ) 
F10.7  ) 


lb  FORMAI  (  1SX »IS»1SX*F10.7*17X  »E15.8  *  15X*  F10.7  ) 

20  CONTINUE 

24  XSEC  =  FOUKPI  *  0(1) 

WRITE  (br 23)  XSEC 

25  FORMAI  (///// 30X2 7H INTEGRA TED  CROSS  SECTION  =  E15.8*  bH  BARNS  //) 
RETURN 

ENU 


ANG00030 
ANG00040 
ANG00050 
ANG00060 
ANG00070 
ANG00080 
ANG00090 
ANG00100 
ANG00110 
ANG00120 
ANG00130 
ANG00140 
ANG00150 
ANG00160 
ANG00170 
ANGOOloO 
ANG00190 
ANG00200 
ANG0021 0 
ANG00220 
ANG00230 
ANG00240 
ANG00250 
ANG00260 
ANG00270 
ANG00280 
ANG00290 
ANG00300 
ANG00310 
ANG00320 
ANG00330 
ANG00340 
ANG00350 
ANG00360 
ANG00370 
ANG003B0 
ANG00390 
ANG00400 
ANG00410 
AUG00420 
ANG00430 
ANG00440 
ANG00450 
ANG004o0 
ANG00470 
ANG004B0 
ANG00490 
ANG00500 
ANG00510 
ANG00S20 
ANG00530 


SUBROUTINE  CftEFS 


Purpose:  Computes  Clebsch-Gordan  coefficients. 

Method:  For  computer  computation  economy,  the  requisite  square  roots  of 

integers,  factorials  and  square  roots  of  factorials  are  stored 
in  common.  The  routine  uses  explicit  algebraic  expressions  for 
Clebsch-Gordan  coefficients  whose  smallest  angular  momentum  value 
is  two  or  less  (with  appropriate  p  tation  of  indices).  Other¬ 
wise  the  general  formula  is  used,  with  a  specialization  for  the 
parity  Clebsch-Gordan  coefficients.  The  input  variables  are 
double  the  angular  momenta  quantities ,  in  order  to  use  them  in 
integer  mode. 

Subroutine  called:  None 


Subroutine  called  by:  FILL,  HUM,  MUSS 
Labelled  c  ommon:  FAC 

Argument  sequence:  (ji,  J2,  J3,  Ml,  M2,  C) 
Argument  List: 


Dimension 


Meaning 


Jl,  J2 ,  J3 


I  Angular  momenta 


Ml,  M2 


I  Magnetic  quantum  numbers 


C 


R 


Clebsch-Gordan  coefficient 


COEFS 
Page  1  of  2 


Enter  CQEFS  | 


Enter  with  J]_,  Jj,  J3, 


=  Mx  + 


Is  min  J  >  2 


Symmetry  operations  to  put  min  J  in  second  place 
with  the  corresponding  M  positive  definite 


Conpute  C  from  appropriate 
special  algebraic  formula  (l) 


Return 


Is  =  0 


Compute  C  from  general  formula  (2) 


Return 


“*** 


Compute  factor  to  reduce  to  parity  C  (3) 


■  Compute  C 

from  formula  for  parity  C  (3) 

References: 

(1)  E.  V.  Condon  and  G.  H.  Shortley,  Theory  of  Atomic  Spectra,  (Cambridge 
University  Press,  1935). 

(2)  G.  Race.,  Phys.  Rev.  438  (1942). 

(3)  L.  C.  Biedenharn  and  M.  E.  Rose,  Rev.  Mod.  Phys.  25>  729  (1953). 


»lBMt  COt> 

SUBROUTINE  COEFS  (  Jl# J2# J3#MI #M2#C  ) 
COMMON  /F AC/F (67 > »KT(9b) #H (bO ) 

M3  s  Ml  ♦  M2 
C  =  U.U 
SI6M  =  l.o 

JMIN  2  M1NU  (J1»J2»J3> 

IF  (  JMIN  *GT*  A  )  60  TO  600 

IF  <  JMIN  .t&.  J2  )  GO  TO  220 

IF  t  JMIN  *EQ.  J3  )  60  TO  230 


210  LI  =  J2 

L2 

=  J1 

L3 

=  J3 

LM1 

=  -M2 

LM2 

=  -Ml 

LM3 

=  -M3 

GO 

10 

220  LI  =  J1 
L2  =  J2 
L3  =  J3 
LM1  =  Ml 
LM2  =  M2 
LM3  -  M3 
60  10  240 

23U  LI  =  Jl 
L2  =  J3 
L3  =  J2 
LM1  =  Ml 
LM2  -  -M3 
LM3  =  -M2 

blGM  =  R IL2+1 ) /K (L3+1 ) 

IF  l  MOU  (  Jl  -  Ml  #  4  )  *NE.  0  )  SIGM  S  -  SIGM 

24U  IF  l  LM2  )  24b#  2b0.  2b0 
24b  LM1  =  -  LM1 
LM2  =  -  LM2 
LM3  =  -  LM3 

IF  i  MOD  (  LI  ♦  L2  -  L3  »  4  )  .NE*  0  )  SIGM  =  -  SIGM 

2bU  JMIN  =  JM1N+1 
K  -  L 1+LM3 
L  =  L1-LM3 

GO  TO  l2bb#26U#30Q#7UU#400> #  JMIN 

2bb  IF  1L1-L3)  B00#2bb#800 
2bt>  IF  (LM1-LM3)  800#2b7#B00 
2b/  C  =  blGM 
GO  10  BUU 

2G0  IF  (L3-L1-LM2)  265»280»270 

2Gb  SIGM  =  -blGM 

2/0  K  =  L 

28U  KF1  =  K  ♦  1 

C  =  K(KHl)/K(2*Ll+2> 

2V0  C  =  blGM*C 

GO  10  800 

3UU  IF  (L3-L1)  30b»310»31b 
30b  IF  «LM21  800»32b#330 
310  IF  (LM2)  800*340# 34b 
31b  IF  (LM2)  80U#3bb»36U 
32b  C  =  — KtL)*KlK)MR«24Ll)*H(Ll  +  l)  ) 

GO  TO  380 

330  C  =  H(L)*K(L+2)/(2.0*H(Ll)*K(Ll+l)> 

GO  10  380 

34U  C  =  FLOAT (LM3)/(H(L1)*K(L1F2) ) 

GO  10  380 

-  U2  - 


COFOOOIO 
COF00020 
C0F00030 
COF00040 
CCF00050 
COF00060 
COF00070 
COF00080 
COF00090 
COFOOIOO 
COFOOllO 
CCF00120 
COF00130 
COF00140 
COF00150 
COF00160 
COF00170 
CO^OOlbO 
COF00190 
COF00200 
COF00210 
COF00220 
CGF00230 
COF00240 
COF00250 
COF002G0 
COF00270 
COF002d0 
COF00290 
COF00300 
COF00310 
COF00320 
COF00330 
COF00340 
COF00350 
COF00360 
COF00370 
COF00360 
COF00390 
COF00400 
COF  00410 
COF00420 
COF00430 
COF00440 
COF  00450 
CGF0P460 
COF00470 
COF00460 
COF0P490 
COF00500 
C0F0P510 
COFOC520 
COF00530 
COF00540 
COFOOSbO 
COF005G0 
COF00570 
COF00580 
COF00590 
COFOOGOO 
COF00610 


€ 


34b  LOX  =  L+2  COF0062O 

t  X  -K(KI«K(L0X)/(H(2*L1>«H(L1*2>>  COF00630 

00  10  380  COF00640 

3bb  LOX  =  L*2  COF00650 

LAX  =  K*2  COF00660 

C  x  N(LOX) *K(LAX) /(R(2*Ll+2)*K(Ll+2) 1  COF00670 

00  10  38U  COF00680 

300  C  =  K(K)*R(K+2>/(2.0*R(Ll+l)«H(Ll«2) )  COF00690 

38U  t  s  biGM«C  COF00700 

00  10  8UU  COF00710 

400  M  =  LM2/2+1  COF00720 

J  =  (L3-L1J/2  +3  COF00730 

00  TO  (480*510  »S40) »  M  COF00740 

480  00  10  (48S*49Ut49b»bOO*SOb) *  J  COF00750 

48b  L  =  R(3)«R(L)*N(L-2)*R(IO*R(K-2)/(R(8>*R(Ll-2)*R(Ll-l)*R(Ll)  COF00760 

1*R(L1+1 ) )  COF00770 

00  1C  b7b  CGF00780 

49U  C  =  -0.b*FL0Al (LM3)*K(6)*R(L)*R(K)/(R(Ll)*R(Ll-2)*R(Ll+l)*R(Ll+2) )COF00790 
bO  10  b7b  COFOO0GO 

4Sb  L  =  O.b  *FLUA1 (3*LM3*LM3«Ll*(Ll*2>>  COFOCeiO 

1  /(R(Ll)*R(Ll-l)*R(Ll+2)*R(Ll+3) 1COF00820 

b7b  COF00830 

COF0084O 

LAX  =  K+2  COFG0850 

C  =  U.b*FL0AT(LM3)  COFOOBOO 

L  *R(6)*R(L0X:-n(LAX)/(R(Ll)*RCLl+l)*R(Ll+2)*R(Ll+4) )  COF0087C 

b75  COF00860 

COF00890 

LAX  =  L+2  COF00900 

LLX  =K+4  CGF00910 

LXX  =  K*2  COF00920 

C  =  H(3)  *H  (LOX )  *H  ILAX )  *HILLX )  *R  (LXX  )/(R(8)*R(Ll+l)  *R (Ll+2)  *R (Ll+3) COF00930 
1*R(L1*4J)  COF0C940 

b7b  COF00950 

(blb*b2U»b2b»530*b3b) *  J  COF00960 

bib  t  s  -K(L+2)*H(L>*R(L-2>*R(K-2J/(2.0*K(L1-2>*R(L1-1)*R(L1>*R(L1'U> JCOF00970 

b7b  COF00980 

(T(Ll+2*LR3-2)  COF00990 

♦R(L*2)*R(L)/(R(LI)*R(L1-2)*R(L1-U)*R(L1*2) >  COFOIOOO 
b7b  COFOIOIO 

b2b  LOX  =  L+2  COF01C20 

L  =  (1.0-FLOATtLM3)l*K(3)*K(LOX)*R(K)/lR(2*Ll)*R(Ll-l)*R(Ll+2)  COF01030 

1*R ( Ll+3) )  COF01040 

b7b  COFOlObO 

i*LM3-Ll-4)  COFOlOoO 

*R(K+2J*R(K>/(R(L1>*R(L1U)*R(L1*2)*R(LU4)  )  COF01070 
b 7b  COFOlOdO 

COF01090 

(.=H(LOX)*R(K+4)*R(K+2)*R(K)/(2.0*R(Ll'U)*R(LU2)*R(Ll*3)*R(Ll+4)  )  COFOUOO 
b  7b  COFOlllO 

(b4b»bbU»bbb*bbU»b6b)  *  J  COF01120 

b4b  L  =  K(L-2)*H(L)*R(L+2)*H(L+4)/(4.0*R(Ll)*R(Ll-2)*R(Ll-l)*K(Ll+l) )  COF0U30 


00  10 

bUU  LOX  =  L+2 


GO  10 

bUb  LOX  =  L*4 


60 

biu 

GO 

bib 

C 

00 

b2U 

C 

1 

GO 

GO 

TO 

b3U 

L=, 

1 

b*FLOAT 

GO 

10 

b3b 

LOX 

=  L+4 

L-R 

(LOX ) *R 

00 

10 

b4U 

00 

10 

b4b 

L  x 

K(L-2> 

GO 

10 

bbU 

L= 

-K(K-2) 

00 

10 

bbb 

LOX 

=  L+2 

LAX 

=  L+4 

L  = 

R(3)*R 

1*H  ( 

Ll+3)  ) 

00 

10 

bbU 

LOX 

-  L+4 

L  = 

-R(K-2 

bib 


COF01140 


bbU  L=  -H(K-2)*K(L)*R(L+2)*K(L+4)/(2#0*R(Ll)*R(Ll-2)*R(Ll+2)*R(Ll  +  l)  )  COFOUbO 

b  7b  COFOUbO 

COFO 1 170 
COFOUOO 

■2>*R(K)*KlL0X)*H(LAX>/(R(8>*R(LD*R(Ll-l>*R(Ll  +  2)  COF0U90 

COFOUOO 

b  7b  COF01210 

COF01220 

L  =  “K(K“2)*K(K)*K(K+2)*R (LOX )/(2»O*R(Ll)*R(Ll+2)*R(Ll+4/*R(Ll+l))COFO1230 
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r  n 


CO  TO  575 

bbb  C=K(K-2)*KU)*R(K+2)*N(K+4)/(4.0*R(Ll  +  l)*R(Ll+2)*R(Ll+3)*R(Ll+%) ) 
5/5  t  =  blGM*C 

GO  TO  800 

/UU  M  =  ILM2+11/2 
J  =  (Lb-Ll+5)/2 
GO  TO  i  710#  740  )  *  M 

710  GO  10  (  72U •  725 »  730*  735  )  »  J 

/2U  t  =  H(3)*KlK-l>*RlL-i)*K(L+l)/(R(8)*R(Ll)*R(Ll-l>*R<Ll+l> ) 

&0  10  780 

/2b  t  =  -FL0ATULl+3+LM3-l)/2)*RtL+l)/(Rt2)*R(Ll-l>*RCLl+l>*R(Ll+2>> 
bO  10  780 
7bU  KP1  =  K  ♦  l 

t  =  -FLOATl 13-3*LM3+L1)/2)*R{KP1)/(R12)*R(L1)*R(L1+1)*R(L1+3)1 
bO  10  780 
7bb  LOX  =  L+b 
KP1  =  K  ♦  1 

KPb  =  K  ♦  b 

t  =  Klb)*KlKPl)*RlKH3)*K(L0X)/(R(8>*RtLl+l)*RCLl+2)*RCLl+3)) 

GO  10  760 

/HU  60  iO  I  7bU*  75b*  760*  76b  )  »  J 

/bu  t  =  -KlL-l)+KlL+l)*K<L+3)/(R<8>*R(Ll>*R(Ll-l>*R(Ll+l) ) 

GO  10  780 
7bb  LOX  =  L+b 
LAX  -  L+l 

C  =  K lb) *R IK-1 )  *R (LAX ) *K (LOX) / (R (8) *R  (Ll**l)  *R  (Ll+1 ) *R (Ll+2) ) 

GO  10  78U 
760  LOX  =  L+3 

t  =  “Nib)  *K IK“1 ) *R (K+l )  *R  (LOX)  /'  (R  (8)  *R  (LI  )*R  (Ll+1 ) *R  CL1+3) ) 

GO  10  780 

7bb  C  =  NlK-l)*H(K+l)*HlK+3)/lK<8)*RlLl+l)*R(Ll+2>*RlLl+3) ) 

/BU  C  =  t+bIGM 
GO  10  800 

THlb  lb  THt  BEGINNING  OF  COMPUTATION  OF  C “COEFFICIENT 
UbING  I  ME  GENERAL  EXPRtbSION. 
bUU  LI  =  J1+J2-J3+1 
L2  ~  J1-J2+J3+1 
Lb  =  -J1+J2+J3+1 
L10  =  J1+J2+J3+3 

IF  l  Mb  .EG.  0  )  GO  TO  615 

BUb  LH  =  Jl+Ml+1 

Lb  =  Jl-Ml+1 

Lb  =  U2+M2+1 

L  7  =  J2-M2+1 

Lb  =  Jb+Mb+1 

L9  =  Jb-Mb+1 

bT  =  KTILIU)  /  (  HTIL1)*RT(L4)*RT(L5>  ) 

bT  =  bT  /  l  RT(Lb)  *  HTIL7)  ) 
bT  =  bT  /  l  R ( J3+1 ) *RT (L2) *RT (L3)  ) 
bT  =  bT  /  l  KT (L8)  *  RT (L9)  ) 

N7  =  L1-L7 
N4  =  Ll-LH 
MIN  =  MAXU  l 0  t N4  »  N /) 

MAX  =  MINU  (L1*L5*L6) 

IF  (M0U(MIN*4) .NE.U)  bIGM=-1.0 

MIN  =  MIN+1 

N1  =  Ll+1 

Nb  =  Lb+1 

N6  =  Lb+1 

bUM  =  U.U 

UO  biu  LZ-NilN* MAX  * 2  .  , 

NIL  =  Nl-LZ  -  '|l!  - 
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C0F01240 

C0F012S0 

COF01260 

C0F01270 

COF01280 

C0F01290 

COF01300 

C0F01310 

COF01320 

C0F01330 

COF01340 

COF013SO 

COF01360 

C0F01370 

COF01380 

C0F01390 

COF01400 

COF01410 

COF01420 

C0F01430 

COF01440 

COF01450 

COF01460 

C0F01470 

COF01480 

COF01490 

COF015C0 

COFOlSlO 

COF01520 

COF01530 

COF01540 

COF01550 

COF01560 

COF01570 

COF01580 

COF01590 

COF01600 

COF01610 

COF01620 

COF01630 

COF01640 

COF01650 

COFOlboO 

COF01670 

COFO160O 

COF01690 

COF01700 

COF01710 

COF01720 

COF01730 

COF01740 

COF01750 

COF01760 

COF01770 

COF01780 

COF01790 

COF01800 

COFOlSlO 

COF01820 

COF01830 

COF01840 

C0FC1850 


NbL  =  Nb-L* 

NbL  =  hfe-LZ 
NbL  =-N4*LZ 
N/L  =“N7*LZ 

T£H*  =  ST  *  F(LZ>  *  F  (NIL)  *  FlNbL)  *  F(N6L)  *  FCNM.)  *  FtN7L> 
C  =  C  ♦  blbM  f  TERM 
610  S16M  =  -blbM 
bO  10  BOU 

bib  IF  t  M2  .tU.  U  )  60  TO  620 

IF  1  lABb(M2)  .NE.  2  )  GO  TO  60b 

JMOU  =  MOO  (  (L10  ♦  1  )»  4  ) 

IF  t  JMOU  .Nt.  0  )  GO  TO  60b 

bT  =  J3  *  IJ3  ♦  2)  -  J1  *  (J1  ♦  2)  -  J2  *  ( J2  ♦  2) 

S1GM  z  U.b  •  blbM  *  ST  /  (  R(Jl)*R(Jl  ♦  25*R( J2)*R(J2  ♦  2)  ) 
620  JMOU  =  MOO  1  (  LI  •  1  )•  fl  ) 

IF  <  JMOU  .Nt.  0  )  SIGH  =  -  SI6M 

L4  =  l  LI  ♦  1  )  /  2 

Lb  =  l  L2  ♦  1  )  /  2 

Lb  =  i  Lb  ♦  1  )  /  2 

L7  =  1  L10  -  1  )  /  2 

C  =  blbM  *  RIJ3  ♦  1)  *  F(L7)  /  (  F(L4)*F(L5)*F(L6)  ) 

L  =  L  *  KIIL1)  *  Kf(L2)  *  KT(LJ>  /  RTCL10) 
ttUU  KtlUKN 
tNU 


COF01860 
COF01870 
COF01880 
COF01890 
COF01900 
COF01910 
COF01920 
COF01930 
COF01940 
COF01950 
COF01960 
COF01970 
COF01980 
COF01990 
COF02000 
CGF02010 
COF 02020 
COF02030 
COF02040 
COF02050 
COF02060 
COF02070 
COF0208G 
COF02090 
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SUBROUTINE  DERIV 


Purpose:  Computes  for  the  Runge-Kut  a  integration  the  derivatives  of 

the  bound  state  vavefunctions  and  the  integrands  of  the  matrix 
elements  up  to  one-half  Bohr  radius. 

Method:  Calculates  the  derivative  of  the  radial  components  from  the 

coupled  Dirac  radial  equations.  In  evaluating  the  in t errand 

y 

of  the  matrix  elements  the  r  factor  is  restored  if  r  <  l. 
Subroutine  called:  SPHBES 
Subroutine  called  by:  RKUT,  RADINT 

Variables  in  unlabelled  Common:  PI,  HALFPI,  F0URFI,  RAD,  SQ£,  Q,  ZA,  ZAZA, 

EFN,  EGN,  V,  CG,  GAM 


Labelled  Common:  BESSEL, 

DF.FUIIC, 

LIMIT,  TAPES,  VECT 

Local  Variables: 

Name  Dimension 

Mode 

Meaning 

Z 

R 

Photon  momentum  *  radius 

i: 

I 

Indexing  variable 

A 

R 

Sum  of  gammas  of  bound  and  free  state 
electron 

Mo 


sibmc  Dfc.ni 

bUBKOUT INt.  ulh:v  deroooio 

L  COMPUTEb  DtHi VAT  IVES  OER00020 

COMMON  Pl#hALFPl»FOURPl#HAO#S«2»Q»2A»2A2A.EFN»EGN»V>C6<30)»GAM<30)DER00030 
COMMUN/btbbtL/FL(lb>»PC(15>»OF(lb»15>»Ml,M2»B(15)  DER00040 

COMMON/OFUNC/F  C  30 ) » 6 ( 30 1 * OF 1 30  > » DG ( 30  > » DFK ( 200 ) #  DFKP ( 200  > » CF  ( 30 ) » HDER00050 
COMMON  /L1MIT/JM»LM»KM#K2M»1END>NE*»NK>NKP»JKB>LMKB>NTAB  DER00060 

COMMON  /TAPtb/X(lbOO)»bCF (1500) »FB<1500>. 6B(1500> »GAMB»SCREEN  DER00070 

COMMON/VtCI /KF (20U) »K6(200 ) »LbES(200) *LBS(200) »LKB  DER00080 

IF  (N1AB«GI*3)  60  TO  10  DER00090 

=  1  DER00100 

00  10  bU  DER00110 

10  IF  (NErf.EQ.O)  60  TO  30  DER00120 

Nt*  =  0  DER00130 

V-“bCF (NT Ad ) /X (NTAB)  DER0014C 

l  -  U*X(NTaB)  DER00150 

CALL  bPHHtb  (2)  DER00160 

00  10  31  DER00170 

3U  NE*  =  1  DER00180 

31  UO  3b  N— 1  # K2M  DER00190 

uhn)  =  cmn)*f(n)/x(ntab)-iefn-v)*G(n)  DER00200 

3b  U6(N)  =  CG(N)*6(N)/X(NTAB>+(EGN-V>*F(N>  DER00210 

UO  HU  N=1#NK  DER00220 

1  =  K61N)  DER00230 

L  =  LbtblNI  DER00240 

UFA ( N ) =b (L ) *6 ( I > *Fb ( NTAb )  DER00250 

IF  (X(NTAB) .6T.1.U)  60  TO  40  DER00260 

A  =  GAMdMGAMB  DER00270 

UFA  IN)  =  UFK(N)*(X(N1AB> )**A  DER00280 

*♦0  CONTINUE  DER00290 

UO  **1  N-l » NKP  DER00300 

1  =  nMNi  DER00310 

L  =  LdblN)  0ER00320 

UFKH(N)=BIL) *F( 1 ) *6b(NTAB)  DER00330 

IF  (X(NTAB) . 6T • 1 *0 )  60  TO  41  DER00340 

A  =  6AM l 1 ) ♦GAMB  DER00350 

UFAPIN)  =  UFKP(N)*X(NTAB)**A  DER00360 

HI  CONIiNUE  DER00370 

bU  KETUKN  DER00380 

EUU  DER00390 
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SUBROUTINE  FILL 


i 


Purpose:  Computes  0(k,X) 

Method:  0(k,x)  is  computed  in  two  passes,  one  for  K  (k  h')  contributions 

the  other  for  K^(k'  h)  and  then  the  two  pieces  are  added  in 
subroutine  HUM  to  yield  0.  The  selection  rules  are  examined 
for  each  matrix  element  contribution  to  determine  what  values 
of  (\  -  t)  are  allowed.  If  X  =  i  ±  1,  f  can  assure  only  the  one 
value  f  =  X  +  l/2 ;  if  X  i  both  the  preceding  values  of  f  are 

possible.  In  the  final  run-through  the  remaining  selection 
rules  required  for  R  are  checked,  R  is  computed  from  explicit 
algebraic  expressions  for  the  Racah  coefficient  with  =  l/2, 
the  corresponding  Clebsch-Gordan  coefficient  called  and  0(h,X) 
formed  with  appropriate  l  sumnation.  The  indices  used  for  0 
are  X  and  K  (K  is  a  positive  integer  uniquely  related  to  h). 
Subroutine  called:  C0EFS 
Subroutine  called  by:  HUM 
Labelled  Common:  FAC,  LIMIT,  QUANT 
Argument  sequence:  (TK,  KW,  LB,  NT,  LP,  Fl) 

Argument  List: 


Name 

Dimension 

Mode 

Meaning 

TK 

200 

R 

Matrix  elements 

KW 

200 

I 

Index  for  k  values  in  matrix  elements 

LB 

200 

I 

Photon  angular  momentum  +  1  for  matrix 
elements 

NT 

I 

Number  of  matrix  elements  K.(h  h’)  or 

K^k’k) 

LP 

I 

/  .  or  l  , 

"K  rV 

FI 

50  x  15 

R 

0+(k,X)  or  0_(h,X) 
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■tt  i.n. 1  fj 


Local  Variables: 


K 


I  Index  of  free  electron  state 


ID  I  Twice  bound  electron  l 

+ 

N  I  Index  of  matrix  element 


L 

JC 

NA,  NB 
LAM 
LEF 
NL 

IR,  IS 
R 


C 


I  Riot on  angular  momentum 

I  Twice  free  electron  j 

K 

I  Loop  index  for  the  range  of  -  l) 

I  2*1 

I  2  *  (\  ±  1/2 ) 

I  (X  -  l)  +  2 

I  Terms  in  R 

R  R  UH,  i-K' 

R  Clebsch-Gordan  coefficient 


JPL  I  2  x  j  ,  +  2  x  j  -2x  photon  angular 

H  H 

momentum 

JML  I  j  2  x  1,  *  2  x  j  j  -  2  x  photon  angular 

momentum 


c 


Enter  TUIt 


E 


D 


Loop  Khrcn&(3) 


f  JSZ  \ 

1  Of  2j 


0 - H  mtua..;±  (SxT^~] 


( 


£ 


Is  photon  4  >  max  4  V 

~~T 


Yes  ^  t  Eliminates  high  4  terns 
in  reduction  loop 


No 


For  each  matrix  element,  examine  the 
selection  rules  [d(4IX),  X^l,  d  (.j^tX)]  toj 
index  the  permissible  values  of  X  -1+2 


© 


Loop  oni-i+2 
tc/2Ti 


Is  X  ■  4  +  1 


j£fi_ 


I 


*© 


No 


f  ■  X  +  1/2 

I 


Compute  R  from  appropriate  algebraic  expression 
(according  as  &«<  4^») 


Multiply  by  matrix  element  and  C-coefficient  and 
add  into  (h,X) 


L 


© 


< 


I 


Is  X  «  4  -  1 


I 


)  **  *® 


No 


f  =  X  -  1/2 


5 
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SlttML  FILL 

bUBKOUUNt  FILL  (TK *KW »LB#NT #LP#Fl ) 

COMMON  /FAC/FACT (67) * RTF AC (96) # ROOT (50) 

COMMON  /LlMlT/JM#LM#RM#K2H#ItNO»N£W#NK#NKP» JKB#LMKB»NTAB 
COMMON/ QU AN I /LF (30) »LMM  30 . JK (30) #FKAP(30 ) #SN(30) #SI (30) #CR(30) 
UIMtNSlON  IK(200)#KW(20U) »LB(200>  »FI(3Q»  15) 

LMF=LMM 
DO  3  LA=1»LMH 
OO  3  K=1.K2M 

3  FIU»LA)=U.U 
LO=2*LF» 

UO  22  N  =  If  NT 
L=LB(N)-1 

IF  (L.GT  »LM)  GO  TO  22 

K=M»(N) 

JC-JK (K ) 

IF  (L.tQ.O  GO  TO  8 

L2=2*L 

JFL=JKB+JC-L2 
JML=IABb  (JKB-JO-L2 
IF  (JFL.LT. 2)  GO  TO  9 

4  IF  (L.LQ.l)  GO  TO  b 

IF  (JML.GT. (-2) )  GO  TO  6 

NA=i 

NB-J 

GO  TO  14 

6  IF  (JML.GT.O  GO  TO  B 

NA=2 
NB=J 

GO  lU  14 
b  NA=J 
N8=J 

GO  10  14 

9  IF  (L *t0. 1 )  GO  TO  13 

IF  (JPL.Gt.U)  GO  TO  11 

NA=1 
NB-1 

GO  10  14 

11  IF  (JML.GT. (-2) )  GO  TO  13 

NA=1 


NB-2 

GO  10  14 

13  NA-2 
NB=2 

14  UO  21  NL=NA»Nb 
LA=L-2+NL 
LAM=2*LA 

IF  (NL.tQ.J)  GO  TO  lb 

LLF-LAM+l 

IF  ( (LU+JC) .LT.LEF)  GO  TO  17 

IF  (lABb(LL)-JC).GT.LtF)  GO  TO  17 

IF  (JKb.GT.LD)  GO  TO  lb 

1K-JKB+ JC+LAM+4 

lb=JKB-JC+LAM+2 

K=-KOOI ( IW) *KOOT ( lb) *KOOT (LA ) 

GO  10  lb 

lb  1S=JC+JKB-LAM 

IF  Ub.tU.U)  GO  TO  17 

1K-JC-JKB+LAM+2 

H=KOOT(IH)*HOOT(Ib)*KOOT(LA) 

16  CALL  COtFb  ( LtF  < LU » JC  » 1 »  U  » C ) 
FllK»LA)=FllK#LA)+TMN)*K*C 
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FIL00010 
FIL00020 
FIL00030 
FIL00040 
FIL00050 
FIL00060 
F1L00070 
FILCOOBO 
FIL0C090 
FIL00100 
FIL00110 
FIL00120 
FIL00130 
FIL00140 
FIL00150 
FIL00160 
FIL00170 
FIL00180 
FIL00190 
FIL00200 
FIL00210 
FIL00220 
FIL00230 
FIL00240 
FIL00250 
FIL00260 
F IL00270 
FIL00280 
FIL00290 
FIL00300 
FIL00310 
FIL00320 
FIL00330 
FIL00340 
FIL00350 
FIL00360 
FIL00370 
FIL00380 
FIL00390 
FIL00400 
FIL00410 
FIL00420 
FIL00430 
FIL00440 
FIL00450 
FIL00460 
FIL00470 
FIL00480 
FIL00490 
FIL00500 
F1L00510 
F1L00520 
FIL00530 
FIL00540 
FIL00550 
F IL00560 
FIL00570 
FIL00580 
FIL00590 
FIL00600 
FIL0Q610 


lr  ^  <NL.£tt,l)  so  TO  21 
Itt  LtF=LAM-l  °  dl 

IF  » tL04JC).LT.LEF)  GO  TO  21 
IF  UABS(LU-JC)  .6T.LLF)  GO  TO  21 


IF  TJKb.GI.LD)  GOTO  19 
lb=JC-JKb*LAM  19 

IF  tlb.Ett.U)  GO  TO  21 
lR=JOJKB-LAM+2 
GO  10  20 

IV  lS=JKb-JG+LAM 

IF  llb.EQ.u)  GO  TO  21 
lK-JKb+JC+LAM+2 

2U  K=NOOT(IH)*ROOT ( lb) *N001 (LA+1) 
J-ALL  COfcFb  ( LEf  » LU » JC  » 1 » (J » C  ) 

FI  (K*LA)=T-1(k»LA)+TK(N)*K*C 

^1  CONTINUE 

n  coniinue 
ketukn 
enu 


60  TO 


GO  TO  21 


F1L00620 

FIL00630 

FIL00640 

F1L00650 

FIL00660 

FIL00670 

F1L00680 

FIL00690 

F1L007Q0 

FIL00710 

FIL00720 

FILCD730 

FIL00740 

FIL00750 

FIL00760 

FIL00770 

FILOO  *0 

FIL0079C 

F IL00800 
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SUBROUTIHE  HUM 


C 


Purpose:  Computes  H(h,u') 

Method:  For  every  h  >  0,  the  selection  rules  on  the  Clebsch-Gordan 

coefficient  we  examined  to  determine  the  allowed  range  of  X 
values.  The  maximum  (positive)  value  of  p*  and  its  minimum 
(largest  negative)  value  are  computed  and  indexed  for  subsequent 
use.  Since  j  depends  only  on  the  absolute  value  of  h}  the 
Clebsch-Gordan  coefficient  is  the  same  for  (-h)  as  for  h.  Over 
the  allowed  range  of  positive  p',  the  Clebsch-Gordan  coefficients 
are  computed  for  given  X  and  h  ,  multiplied  in  turn  by 
0(h,\)  and0(-K,k)  to  obtain  corresponding  contributions  to 
H(h,p’)  and  H(-x,p'),  and  the  products  are  summed  over  X.  The 
results  are  denoted  by  HF(K,M)  where  K  is  a  positive  integer, 
indexing  h,  and  M  =  p '  +  l/2.  The  procedure  is  then  repeated 
for  p'  negative,  leading  to  HFM(K,M),  where  K  is  as  above  and 
M  =-p'  +  l/2.  The  explicit  separation  of  positive  and  negative 
p'  terms  is  useful  later  on. 

Subroutines  called:  C0EFS,  FILL 

Subroutine  called  by:  RABINT 

Labelled  Common:  FID0,  LIMIT,  MAT,  QUANT,  TRANS,  VECT 

Local  Variables: 

Name _ Dimension _ Mode _ Meaning _ 


L 


IMP 

K 

LA 

JC 

LAD 


I  Maximum  photon  angular  momentum  +  1 

I  Loop  index  for  bound  electron  state 

X 

I  Twice  j  (free  electron  state) 

I  Minimum  X  from  selection  rules 
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Name 


Dimension 


Mode 


Meaning 


LAP 

JPJS 

JKEG 

M 

LAM 

MU 


I  Maximum  \  from  selection  rules 

I  * )jnnjc  +  provided  u'  >  0  permitted, 

-1  otnervri.se 

I  (-u’)  +  1/2 

v  'max  ' 

I  Loop  index  for  p. ' 

I  Twice  X 

I  Twice  u' 


-  to  - 


SlBMt  HUMM 

bUBKOUl iNt  HUM 

COMMON/F 1U0/K I ( 30 » 15 ) # D l 30 ) # JMP i NAME , SHELL * QV . c3 » IZ 
COMMON  /L1MIT/JM#LM'*M'K2M'IEN0»N£M»NK»NXP'JKB'LMKB'NTAB 
COMMON/ MAlYbF (30 ) ibG( 30 ) i FK (200 ) i FKP( 200 ) »SFK (200) i SFKP(200 ) »RCUT 
COMMON/QUAN I /LK ( 30 ) »LMK (30) » JK (30) #FKAP(30 ) »SN(30) »SI (30 ) #CR (30 ) 
COMMON/ TKANb/HF  (30* lb)  *HFM(30fl5)  » JN6C30) » JPS<30) 

COMMON/ VtC l/KF(20U) * KG (200 ) * LbEb ( 200 ) i LBS (200 ) *  LKB 
UlMtNblON  Fl(30f  li>)  *FTP(3U#lb) 

LMP=LM-»1 

CALL  FILL  (FKfKGfLBtb»NKfLMKBf FT ) 

CALL  FILL  ( F  KP  f  KF  f  Lbb  f NKP » LKB  f  FTP ) 

UO  11  K=lfK2M 
00  //  LA=1#LMH 

It  FI(K»LA)=F1 (K#LA)FFTP(K#LA) 

IF  (  JM  .tO.  0  )  RETURN 

UO  at  KSlfK2Mf2 
KF  =  K  ♦  1 
JC=JK(K) 

LAO-1 

IF  (JC.tO.JKb)  GO  TO  78 

LAO=lAbb  (JC-JKb)/2 
fti  LAP  =  ( JC+JKB)  /2 

IF  (  LAP  .01.  LMP  )  LAP  =  LMP 

JPObsjC-2 

IF  (JC.LT.4)  60  TO  70 

JPOb-MINO  (JKb.JPOb) 

JPOb- ( JPOb* 1 )  /2 
t 9  JNtfa=MlNU  (JKBf  (JC*2)  ) 

JNtGS ( JNbG+1 )  /2 
JNO(K ) -JNbG 
JPb ( K  J -JPOb 
JNG(KF)SJNtG 
JPb (KP)  = JPOb 
IF  (JPOb. LI. U) 

DO  ttU  M= It JPOb 
HF(K.M)  S  U.O 
8U  MF(KPfM)  =  U.U 
ttl  UO  82  M— 1 f  JNto 
KFM(K.M)  =  0.0 

82  MFM(KP.M)  s  U.O 
IF  (LAU.GI.LAP) 

UO  bb  LA=LAU.LAP 
LAM=2*LA 
IF  (JPOb. U.O) 

UO  83  M= It JPOb 
MU=2*M-1 

CALL  COtFb  (LAMf JKbf JCf2fMUfC) 

HF(KfM)  s  MF(KfM)  +  FI(K»LA)  *  C 

83  HF(KPfM)  =  HF(KPfM)  ♦  Fl(KPfLA)  *  C 

84  UO  bb  MSl.jNtO 
MU=1“2*M 

CALL  COEFb  (LAMf JKbf JCf2»MUfC) 

MFM(KfM)  S  HFM(KfM)  ♦  F 1 (Kf LA)  *  C 
bb  MFM(KPfM)  =  HFM(KPfM)  ♦  FI(KPrLA)  * 
bb  CONI INUt 
at  CONI INUt 
HtTUHN 
LNO 
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00  10  81 


GO  10  87 


00  TO  84 


HUM00010 
HUM00020 
HUM00030 
HUM00040 
HUM00050 
HUM00060 
HUM00070 
HUM00080 
HUM00090 
HUM00100 
HUM00110 
HUM00120 
HUM00130 
HUM00140 
HUM00145 
HUM00150 
HUM00160 
HUM00170 
HUM00180 
HUM00190 
HUM00200 
HUM00210 
HUM00220 
HUM00230 
HUM00240 
HUM00250 
HUM00260 
HUM00270 
HUM00280 
HUM00290 
HUM00300 
HUM00310 
HUM00320 
HUM00330 
HUM00340 
HUM00350 
HUM00360 
HUM00370 
HUM00380 
HUM00390 
HUM00400 
HUM00410 
HUM00420 
HUM00430 
HUM 00 440 
HUM00450 
HUM00460 
HUM00470 
HUM00480 
HUM00490 
HUM00500 
HUM00510 
HUM00520 
HUM00530 
HUM00540 
HUM00550 
HUM00560 
HUM00570 


This  page  is  intentionally  left  blank. 
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SUBROUTINE  INTERP 


Purpose:  Interpolates  on  the  bound-state  tabulations  for  radii  greater 

than  one -half  Bohr  unit  where  the  integration  grid  is  much 
finer  than  the  table,  to  obtain  intermediate  values  of  the  bound- 
state  wavefunctions  and  corresponding  potential. 

Method:  Linear  interpolation  between  successive  entries  in  table. 

Subroutine  called:  None 


Subroutine  called  by:  XDERIV 
Labelled  Comm or.:  LIMIT,  tfNWARD,  TAPES 
Local  Variables : 

Name _ Dimension _ Mode _ Meaning 


DX  R 

DL  R 

QU0T  R 

DV  R 

DGB  R 

DFB  R 


Difference  between  two  successive 
tabulated  radial  values 

Radial  increment;  difference  between 
integrating  radius  and  lower  grid  point 

Fractional  radial  increment 

Difference  between  two  successive 
tabulated  potential  values 

Difference  between  two  successive 
tabulated  values  of  "large"  component 
of  wave function. 

Difference  between  two  successive 
tabulated  values  of  "small"  component 
of  wave function. 
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©■ 


( 


(  Enter  IH7ERP  j 

-^T 


Is  radial  variable  >  ^ _ Yes 

current  radial  grid  value 


[ 


> 


No 


Proceed  to  perform  inte jrpolation 

r 

Compute  difference  of  wo  successive 
grid  radial  values  which  bracket  desirec 
radial  value 

f 

Calculate  radial  increment  and 
fractional  radial  increment 

' 

(potential  values,  major  components  of  wav 
functions  and  minor  components 
of  wave  functions 


Calculate  by  interpolation  new 
potential,  major  &  minor  components 
of  wave  function  corresponding 

_ to.  ra<4ai,  va.1us _ 


15 


Return 
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Sib**  It  1N1H 

SUbMOUUNfc.  iNTtKP  INT00010 

COMMON/L 1 M  i  T /  JM » LM » KM » KiiM »  I ENO » Nfci  »  NK  » MKP »  JKB  *  LMKB » NT A8  INT00020 

COMMON  /ONN AKU/KX »  bC  X » 6b  X # FBX  INT00030 

COMMON/ TAFtb/X l lbUU) «SCF(1500) ’FB(ISOO) *68(1500) *6AMB» SCREEN  INT00040 

b  IF  (KX.Gr.X(NTAbU))  SO  TO  10  INT00050 

OX  s  XlNTAb+1 )~X(NTAb)  INT00060 

OL  =  KX-X(NlAd)  INT00070 

0001  =  OL/UX  INT00080 

OV  =  SCFlNIAb«*l)-5CMNTAb)  INT00090 

06b  =  GB(NTAb+l)*6b(NTAb)  INTOOlOO 

OFb  =  FblNI Ab+1 ) -FbiNTAb)  INT00110 

bCX  s  bCF  l  NT  Ab )  ♦QUOT *l)V  INT00120 

GdX  =  GdlNlAb)+GUOT*UGb  INT00130 

FbX  =  F  b ( N  T  AU ) ♦QUOT  *OFb  INT00140 

GO  10  lb  INT00150 

1U  NTAb  =  NTAo+1  INT00160 

GO  10  b  INT00170 

ib  KEIUKN  INT00180 

fc.NL)  INT00190 


This  page  is  intentionally  left  blank. 
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SUBROUTINE  LEGEND 


Purpose:  Computes  the  Legendre  coefficients  of  the  cross  section. 

Method:  For  j  =  0,  this  consists  of  carrying  out  the  sum 

2  2 

D0  =  ILfi2  £  2\+  1  1  over  all  K  ^  X 

X  K 

(The  0's  have  been  initialized.)  For  j  >  0,  subroutine  MUSS 

is  called  where  the  T.(x.,k’')  of  the  sum 

J 

TT2  2  i 

D.  =  — £  cos (6  -  6  ,,)  T.(K>K")  is  performed. 

0  <*  H  K  j 

The  simulation  (  (-l)J  2  cos(  )r.  (  )  )  is  subject  to  the 

xx"  J 

selection  rules  A(jH<^.tj),  AUH^H.J)  and  iH  +  ZH, ,  +  j  =  even 
integer.  The  diagonal  terms  in  the  double  sum  are  done  first 
(with  the  cosine  equal  to  unity),  carried  over  j  >  1/2  and 
contributing  to  even  j  terms  only.  Since  the  off-diagonal  terms 
are  symmetric  in  x  and  x",  twice  the  sum  with  K"  <  K  is  taken. 
For  given  x  and  x",  the  smallest  j  for  which  there  can  be  a 
contribution  is  I  L  -  l  , ,  provided  I  j  -  j  ,  ,1  is  not  larger; 

*  H  H  "  *  H  H  | 

otherwise  j  =  I  j  -  j  ,,  has  the  wrong  parity  and  the  minimum 
j  value  is  |j^  -  jk,  ,|  +  1  =  1 -  lRl  ,|  +2.  There  may  also  be 
contributions  for  larger  j  (going  up  in  steps  of  two  to  preserve 
parity)  up  to  the  lesser  of  (i  +  i  ,,)  and  (j  +  i  ,,),  or  up  to 

H  H  H  H 

an  assigned  maximum  j  if  smaller. 

Subroutine  called:  MUSS 
Subroutine  called  by:  PELEC 
Labelled  Common:  F3D0,  LIMIT,  QUANT 
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Local  Variables: 


Name 

LMP 

J 

LA 

STAT 

SKP 

JMAX 

JMIN 

K 

TJ 

JDU 

C0D 

MINI 


Mode 

I 

I 

I 

R 

R 

I 

I 

I 

R 

I 

R 

I 


Dimension 


Meaning _ 

Max  photon  angular  momentum  +  1 
Order  of  Legendre  coefficient 
k 

2k  +  1 
Z  02(k ,k) 

K 

Max  j  contribution  for  given  h  and  h 1 ’ 
Min  ^  contribution  for  given  h  and  h" 
Loop  index  for  bound  electron  state 


T  {k,h") 
J 


2  cos  (6k  -  6k,,) 
(-1)0'1 
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9 


(^^7) 


(  18  "  >  nin  j 

- No  ±-Z- 


Set  min  j  *  win  i  +2 

- - I - “ - 


f  Is  mir,  j  >  max  order  of  \ 
V  Leg,  coef.  / 


Yes 


No 


Max  j  s  j  +A,  .  1 1  l  +1  , , ,  JM;increment  min 

H.  H  * 

j  and  max  j  by  1  for  indexing 


Yes 


Is  max  j  <  min  j 


C 


No 


Is  K.E.  for  electron  zero 


Yes 


3 


No 


Compute  2  cos  (6  -6  ,,) 

K  K 


(1) 


■0 


<5 


■0 


Loop  from 
J 


max  toi^l^] 


j  min  to 


© 
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LEG0001C 
LEG00020 
LEG00030 
LEG00040 
LEG00G50 
LEG00060 
LEG00070 
LEG00080 
LEG00090 
LEG00100 
LEG30L10 
LEG00120 
LEG00130 
LEG00140 
LEG00150 
LEG00160 
LEG00170 
LEG00160 
LEG00190 
LEG00200 
LEG00210 
LEG00220 
LEG00230 
LEG00240 
LEG00250 
LEG00260 
LEG00270 
LEG00260 
LEG00290 
LEG00300 
LEG00310 
LEG00320 
LEG00330 
LEG00340 
LEG00350 
LEG003&0 
LEG00370 
LEG00380 
LEG00390 
LLG00400 
LEG0041 0 
LEG00420 
LEG00430 
LEG0044C 
LEG00450 


c 
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AS  MMf 


tlHMC  Ltbt 

bUDKOUTl NE  LEGEND 

COMNON/F IUO/F I ( 30 # 15 ) # D ( 30 ) # JMP # NAME » SMELL # QV # EB # 12 
COMMON/LIMIT/ JM»LM»KM>K2M*1EN0»NCM»NK»NKP»JKB»LMKB»NTAB 
COtiMON/GUANf /LK (39# *LMK (30) * JK(30) »FKAP(30) * SN(30) »SI (30) *CR(30) 
L*P=LMU 
DO  /9  U=1#J*P 

/y  ucjisu.o 

00  01  LASl'LMt* 
bTAT&2*LA+l 
bKK=U.O 
00  80  K=i»K2M 

OU  bKP  =  bKF+F I (K»LA) *F 1 (K»LA) 

01  0  ( 1 )  :L)  ( l )  ♦bKP/bT  AT 

IF  (JM.tQ.U)  GO  10  90 

IF  (JM.EQ.l)  GO  TO  84 

UO  03  K=I»K2M 

IF  (  JMK)  .tU.  i  I  GO  TO  83 

OMAA=MlNU(OK(K) »2*LK(K>  * JM)+1 
DO  02  J=3»JMAX»2 
CALL  MUbb (A»K»J»TJ) 

02  D(J)=D(J)+TJ 
83  C0N11NUE 
04  DO  89  A— 2  » K2M 
KLb=K-l 

DO  80  KK-1 »KLb 
JMIH=lABb(LMK)-LK(KK)  ) 

J01F=lABb( JMK)-JK(KK) ) 

IF  (JDlF.bl .(2*JMIN) )  JMIN=JMIN+2 

IF  ( JMlN.Gl » JM)  GO  TO  38 

JM1N— JMlN+1 

JMAX=M1NUU  (JMK)+JK(KK))/2I »  (LA  (K )  +LK  (KK  ) )  » JM)  +1 

IF  t  JMAX  .LT.  JMIN  )  GO  TO  88 

IF  (  IthO  .tQ.  1  )  GO  TO  83 

COD  =  2.U  *  (  CK (K) *CK (KK)  ♦  SI(K)*SI(KK>  ) 

GO  10  06 
Ob  COD  =  2.U 
06  MINIS  MOD ( JMIN »  2 ) 

IF  (MINI.EU.U)  COOS-COO 
DO  0/  J=JMIN*JMAX»2 
CALL  MUbb  (K  *KK  » J»  1 J) 

0/  D( J)=Ul J)  +  l J*COD 
00  CONUNUL 
09  CONI INUt 
9U  KETUHN 
tND 
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SUBROUTINE  L$GGAM 


Purpose:  Computes  the  natural  logarithm  of  the  gamma  function  for  complex 

arguments,  i.e.  Re  in  T(x  +  iy)  and  I31  in  r(x  +  iy) 

Method:  a)  Set  T(z)  =  -z  +  (z  -  l/2)  inz  +  jtnV&T  +  J(z)  where  j(z)  is 

given  as  a  continuous  fraction.  See  Wall,  "Analytic  Theory  of 
Continued  Functions,"  p.  36U,  formula  93. 9- 
b;  For  x  >  2,  jtn  T(z)  is  computed  from  the  recursion  relation: 
inr(z)  =  jtnT(l  +  z)  -  Inz 

c)  For  negative  x,  the  Tm  jfcnr(z)  can  be  thought  of  as  being  equal 
to  V  +  2nk  where  k  is  an  integer  and  V  is  given  by  this  routine. 
Restrictions:  a)  x  and  y  may  not  both  be  equal  to  zero  . 

b)  if  y  =  0,  x  may  not  be  equal  to  a  negative  integer 
This  routine  is  taken  from  M.  S.  Shapiro  and  M.  Goldstein,  "A 
Collection  of  Mathematical  Computer  Routines,"  NY0-1U80-1U  (1965). 

Subroutine  called:  None 
Subroutine  called  by:  RADINT 
Argument  sequence:  (X,  Y,  U,  V) 

Argument  List: 


Dimension 


Mode  Meanir 


Real  part  of  argument 
Imaginary  part  of  argument 
Real  part  of  result 
Imaginary  part  of  result 


[ 


SlbMC  LGAM 

SUbHOUl iNt  L0G6AM ( X *  Y  *U*  V ) 

C  This  SUbHOUTlNt  COMPUTES  THE  NATURAL  LOG  OF  THE  GAMMA  FUNCTION  FOR 
C  COMPLEX  ARGUMENTS.  THE  ROUTINE  IS  ENTERED  BY  THE  STATEMENT 

*  C  CALL  LOGGAM (XtYtU’V) 

C  KhEKE  X  IS  THt  REAL  PART  OF  THE  ARGUMENT 

C  Y  IS  THt  IMAGINARY  PART  OF  THE  ARGUMENT 

C  U  IS  THE  REAL  PART  OF  THE  RESULT 

C  V  IS  fht  IMAGINARY  PART  OF  THE  RESULT 

DIMENSION  H(  7) 

H(l)=2.2b94«a974 

H(2)=l.bl74/3649 

Hl3)=1.011S230b8 

H(4)=.S25bUb4b9U 

mb)  =  .2S24809S24 

h(b)=0. 0444444343 

H(  71  =  0. 0844344433 

E2=l.b70/9b42b79 

E8=4.141S92bb4S9 

bl=U.U 

1  b2=U.0 

J=2 
X2  =  X 

4  IF  ( X)  1 » 2 * 4 

4  bb=ATAN  (Y/X) 

I  =X*X 

5  b7=Y*Y+T 

t  HEAL  PAKT  OF  LOG 

T 1=0 • S* ALOG ( b7 ) 

IF(X-2.0) 7  *  7 » b 

7  bl=bl+bb 
b2=b2+Tl 
XSX+1.0 
J=1 

GO  10  4 

b  T4=-Y*bb-MTl*lX-.S»-X+.91894bb332  ) 

12=bb*lX-.b)+Y*Tl-Y 
T  4=X 
Tb=-Y 
Tl=b  7 

UO  b  1=1*7 
T=H(i)/Tl 
I  4= 1  * l 4*X 
Tb=“( I ♦Tb+Y ) 

b  11=14*1 4+’Tb*Tb 
I 4=l4«X+  T  4 

j  12=-Ib-Y+T2 

go  io  (y*joi*j 

9  1 4=1 4-b2 

T2=l2-bl 

1U  IF ( X2 ) 1 1  *  12  *  12 
12  U=  T  4 
V=T2 
X=X2 
HE1UHN 
11  U=T4-E4 
V=l2-tb 
X=X2 
RETURN 

C  X  lb  /.LHO 

2  1=0.0 

lFtY)14*l4*lb  vh  _ 


LGM00010 

LGM00020 

LGMC0030 

LGM00040 

LGMOOObO 

LGM00060 

LGM00070 

LGM00060 

LGM30090 

LGMllOlOO 

LGM00110 

LGM0G120 

LGM00130 

LGM00140 

LGM00150 

LGM00160 

LGM00170 

LGM00180 

LGM00190 

LGM00200 

LGM00210 

LGM00220 

LGM00230 

LGM00240 

LGM00250 

LGM00260 

LGM00270 

LGM00280 

LGM00290 

LGM00300 

LGM00310 

LGM00320 

LGM00330 

LGM00340 

LGM00350 

LGM00360 

LGM00370 

LGM00380 

LGM00390 

LGM004U0 

LGM00410 

i.  M00420 

LGM00430 

LGM00440 

LGM00450 

LGM00460 

LGM00470 

LGM004B0 

LGM00490 

LGV00500 

LGM00510 

LGM0P52C 

LGM0053C 

LGM0054C 

LGM0055C 

LGM005bC 

LGM0057C 

LGM005b( 

LGM0059( 

LGM0060I 

LGMOObll 


14  Bb=“t2 
60  TO  b 
lb  Bb=E2 
60  10  b 
X  lb  NtbATlVt 
1  t4=U.U 
tb=U.O 
lEb=U 

lb  t4=t4*.b*(AL0G(X*X  ♦Y*Y  )) 
tb=tb*ATAN  (Y/X) 
lEb=ltb+l 
X=X+1.0 
IF  IX) lb* 17# 17 

1/  IF  C  MOD  ( ltb*2) ) 18*4*16 
lb  tb=tb+t8 
60  10  4 

14  teKl 1  Lib*  IV)  X2»Y 

IV  FOHMAT (2VH  ATTtMPTEO  TO  TAKE  LOGGAM 
CALL  LX  1 1 
HETOHN 
t  NO 


LGM00620 

LGH00630 

LGM00640 

L6M00650 

LGM00660 

LGM00670 

LGM00660 

LGM0G690 

LGM00700 

LGH00710 

LGM00720 

LGM00730 

LGHr'0740 

LGM00750 

LGM00760 

LGM00770 

LGMOO70O 

2HX=F6.0»1X2HY=F6.0)  LGM00790 

LGMOO0OO 

LGM00610 

LGMOO02O 
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SUBROUTINE  MUSS 


Purpose:  Performs  the  p’  sum  and  computes  the  term  T.  (h,h")  = 

1/2,-  1/2)  I  (-l)u*  +  X/2  H(k  ,p  ' )  H(h  1 '  ,p' ) 
p 

C(j  ,,j  j;  p'  +  1,  -p*  -l)  which  is  used  in  computation  of 
H  H 

D  (Legendre  coefficient). 

J 

Method:  The  sum  is  carried  out  first  for  (p'  +  l/2)  ^  0,  computing  the 

C-coefficients  and  recording  them  as  an  indexed  variable,  the 
upper  limit  on  p’  being  the  lesser  of  the  upper  limits  for 
h  and  h1'  computed  in  subroutine  HUM.  The  sum  for  the  negative 
(p*  +  1/2)  is  similarly  carried  up  to  the  lesser  of  the  maximum 
|  p'  J  values  for  negative  p',  but  the  C-coefficient  is  obtained 
from  the  previously  obtained  ones  by  a  symmetry  operation 
(introducing  at  most  a  sign  change).  The  C-coefficient  in  front 
of  the  sum  is  the  one  computed  for  p’  =  -l/2. 

Subroutine  called:  C0EFS 
Subroutine  called  by:  LEGEND 
Labelled  Common:  QUANT,  TRANS 
Argument  sequence:  (K,  KK,  J,  TJ) 

Argument  List: 


KK 


I  Index  for  h ' ' 


J 

TJ 


I  Order  of  Legendre  coefficient  +  1 

R  T.  (h,h") 
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JD 

I 

Twice  order  of  Legendre  coefficient 

C 

R 

Clebsch-Gordan  c  ’efficient 

CL 

30 

R 

Clebsch-Gordan  c oefficient  saved 

JP0S 

I 

Range  of  positive  u'  (-1  means  none) 

M 

I 

Loop  index  for  p' 

PM 

R 

(-1)^  +  1^2  (times  i  K  K  if  n'  negative) 

MP 

I 

2  |  n'  +  1  1 

JNEG 

I 

Range  of  negative  p' 

JM0D 

I 

JK+jK,,+j 

i 
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SlttMC  MUSt 

SUBKOUIINE  MUSS  (KtKK'J'TJ) 

COMMON/ QUANT /LM  30 ) »LMK (30) » JK(30) »FKAP(30) tSNOO) #SI (30) »CR(30) 
COMMON/ TRANS/HF  C  30 ♦ IS ) # HF M ( 30  # 15  > » JN6  C  30 ) » JPS ( 30 ) 

DIMENSION  CL (30) 

J0=2*IJ-1) 

CALL  COtF  S(JK(K) » JK ( KK ) »JD»1»-I*C) 

1J  =  HF  M  l  K  #  1 )  *  HFM(KK’l)  *  C 
CH1)=C 

JPOS=M1NOI JPS(K) # JPS(KK) ) 

IF  IJPOS.tU. l-l) )  60  TO  33 

PM=1.U 

DO  30  M=l#dPOS 

PM=-PM 

MP=2*M+1 

CALL  COtFbl  JMK)  #  JK(KK)  * JO» MP # -MP » C ) 

CL  IMP ) — C 

3U  Id  =  TJ  ♦  hHK*M)  *  HF  IKK ♦  M)  *  C  *  PM 
33  JNtto=MlNUlJNGlK) *JN6lKK) ) 

IF  IJNtG.tO.l)  60  TO  32 

JMOU-MOUt (JK (K>+JKlKK)+JO) »4) 

PM— “1 • 0 

IF  IJMOU.tU.U)  PM-1 • 0 

UO  31  M— 2 » JNtG 
PM— “PM 
MP=2*M-3 

31  1J  =  Id  ♦  tiFM IK » M)  *  HFMiKK  #M)  *  CLCMP)  *  PM 

32  I J= ( J*CL ill 
KLTUKN 

tNU 


MUS00010 

MUS00020 

MUS00030 

MUS00040 

MUS000S0 

MUS00060 

MUS00070 

MUS00080 

MUS00090 

MUS00100 

MUS00110 

MUS00120 

MUS00130 

MUS00140 

MUS001S0 

MUS00160 

MUS00170 

MUS00180 

MUS00190 

MUS00200 

MUS00210 

MUS00220 

MUS00230 

MUS00240 

MUS00250 

MUS00260 

MUS00270 

MUS00280 

MUS00290 


SUBROUTINE  RADINT 


Purpose:  This  is  the  control  subroutine  for  the  integrations. 

Method:  Hie  requisite  coefficients  and  the  initial  values  of  the  free- 

electron  wavefunctions  and  their  derivatives  are  computed.  The 
radial  integrals  are  performed  by  calling  the  Runge-Kutta 
integration  subroutines  in  D0  loops.  Normalization  factors  are 
obtained  frcm  WNORM  and  applied  to  the  matrix  elements.  The 
phase  shifts  are  obtained  by  wave -matching .  Finally,  subroutine 
HUM  is  called  to  start  the  angular  momentum  sums. 

Subroutines  called:  L0GGAM.  RKUT,  DERIV,  XDERIV,  XRKUT,  WN0KM,  SPHBES,  HUM 

Subroutine  called  by:  PELEC 

Variables  in  unlabelled  Gbmmon:  PI,  HALFPI,  ESURPI,  RAD,  SQP,  Q,  ZA,  ZAZA, 

EFN,  EGN,  V,  CG,  GAM 

Labelled  Common:  BESSEL,  DFUNC,  LIMIT,  MAT,  QUANT,  VECT,  INWARD 


Local  Variables: 

Name 

Dimension 

Mode 

Meaning 

E 

R 

2 

Free  electron  energy  (in  me  units) 

RK 

R 

yl?  -  1,  free  electron  momentum 

I 

I 

Loop  index  of  free  electron  state 

TUGAM 

30 

R 

2  s!/  -  ,  Z  ^  +  1 

'137.0367' 

I  UP 

I 

Max  number  of  matrix  elements 

NJS 

I 

Loop  index  to  initialize  to  zero  components 
of  free  electron  wavefunction,  their 
derivatives,  integrand  for  matrix  elements 
and  their  integration  storage  variables. 

N 

I 

Range  on  number  of  iterations  for  radial 
integration  up  to  radius  of  one. 

XCUT' 

R 

Radial  cutoff  on  uniform  grid  in 
tabulation. 

WAVE 

R 

Max  of  wave  numbers  of  free  electron  and 
photon . 
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Name 


Dimension 


Mode 


Meaning 


LUB 


RX 

LTOT 


I  Range  on  number  of  iterations  beyond 

one-half  Bohr  radius  for  radial 
integration 

R  Current  radial  variable  beyond  one-half 

Bohr  radius 
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<£> 


Calculate  no.  of  iterations  on 
Icon  to  integrate  up  to  radius 
- 1 _ r-u  one _ 


Call  RKUT  for  integration 
up  to  radius  =  1 


E 


Insert  remaining  factor  in  initial| 
normalization  of  wave  functions, 
their  derivatives  &  integration 
storage  locations 


Insert  factor  in  matrix  elements 
K  (hx'),  she  derivatives  and 
1  storage  locations 


Insert  factor  in  matrix  elements 
K.(h'k),  the  derivatives  & 

1  storage  locations 


G 


RADINT 

Page  3  of  'bj 


Loop  on  no.  of 
iterations  to 
integrate  up  to.. — . 
radius  of  1  to (350) 

Loop  on  no.  waue 
functions  to 


Loop  on  no. 
matrix  elements 

0K  (kh1)  to 

1 

Loop  on  no.  matrix 
elements 

"K.(h  k)  to  510) 

1  — s 


Ret  step  size  for  integration  beyond| 
radius  of  one,  compute  radial 
cutoff  on  uniform  grid  tabulation 


0 


Call  subroutine  DERIV  to 
supply  derivatives  for  integrator 


I 


Calculate  no.  of  iterations  on  loop 
necessary  to  integrate  up  to 

_ radial  cutafi _ 


- Call  RKUT 

► 

f  Is  cutoff  radius 
l  radial  cutoff 

for  integration 
on  uniform  grid 

D 

,  No 

Must  continue  integration  by 
interpolation  process 

Loop  on  no.  iterations 
to  integrate  uk — \  to 


cutoff  radius  to  1570 


Yes 


Integration 

completed 


3 
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Multiply  matrix  elements  by 
normalization  factor 


Loop  on  no.  matrix 
elements 


Multiply  matrix  elements 

K.(h'k)  by  (-l)(l-K’+Vi)/2 

X 


K  (h'h)  to (401) 

X 


Multiply  matrix  elements  by 
normalization  factor 


Call  subroutine  HUM  to  start 
summing  process 


SltiMt  RAUN 

bUBKOUT iNt  KAUINT  RAD00010 

COMMON  P I  * mAlFPI » FOUHP I » RAO *SQ2*Q»2A» ZAZA # EFN * E6N *  V » CG ( 30 ) » GAM ( 30 ) RAD00020 
COMMON/«tSbtL/FL(lb) »PC(1S) .OF  <lb»15) *M1*M2*B(15)  RA000030 

COMMON/UFUNC/F (30) »G(3U) #DF(30) *06(30) »UFK(200) *OFKP(200) »CF(30) .HRAD00040 
COMMON/LIMIT  /  JM  *  LM  *  KM  *  K2M  #  IENU  *  NEW  *  NK  *  NKP  *  JKB » LMKB » NTAB  RA000050 

COMMON/MAT/bF (30) *bb(30) *FK (200) *FKP(200) *SFK(200) *SFKP(200>  »RCUT  RA000060 
COMMON  /ONNAKO/HX*bCX*6BX*FBX  RA000100 

COMMON /QUANT /LA ( 30 ) *LMK ( 30 ) *  JK ( 30) *FKAP( 30) *SN(30) *SI (30) *CR(30)  RA000070 

COMMON/VtCT/KF (200) »KG(20Q) »LBES(200) »LBS(200) *LKB  RAOOOOBO 

UlMtNblON  HNOKM(30) *  TUGAM(30)  RAD00090 

1UU  FOHMA1 (bX#bOHLtNGTH  UNITS  ARE  HBAR  MC  (  1  BOHR  RADIUS  =  137  )//)RAD00110 
101  FOHMAT  ( SX » 24H INTEGRATION  STEP  SIZE  IS*F11.7»  5X*5HUP  T0»F9.3'/)  RA000120 


t  =  tFN+l.U 
tGN  =  t+l.U 

IF  (  tFN  .to.  0.)  GO  TO  302 
NK  =  SUN  T (LFN*EGN) 

1ENU  =  0 

bdtG  =  bQHT  (tGN) 
bQt  =  KK/bUtG 
NKt  =  HK  /  tGN 
GNU  =  2A*E/NK 
INK  =  2.U*KK 

TwOKK  =  t.XP  (  GNU  *  HALFPI  )  /  SORT  (  TRK  •  PI  ) 
UO  300  1=1*K2M 

GAM ( 1 )  =  buRT  (FKAP( 1 )**2-ZAZA) 
lUGAM(l)  =  2.0*GAM( 1 )+l.U 
CALL  LOGGAM  ( bAM ( 1 ) * GNU  *  XRt  *  X IM ) 

ZK  =  1  WORK* (TRK**GAM(I) )*EXP( XRE ) 

IF  (FKAP(I) .GT.U.O)  GO  TO  268 

CF  ( 1 )  =  FKAPU)-GAM(I) 

CG ( I )  s  -ZAZA/CF(I) 

MM  s  ZAZA  *  EFN  /  (  CF ( 1 )  *  (  E*  FKAP(I)  -  GAM( I ) 
HP  =  2.0  -  HM 
HM  =  bQRT  (  HM  ) 

HP  =  SORT  (  HP  ) 

G ( 1 )  =  ZK  *  SQEG  *  (  GAM ( 1 )  *  HP  +  GNU  *  HM  ) 

F ( I )  =  ZA*G ( I ) /CF ( 1 ) 

GO  10  270 


2btt  (.6(0  =  -FKAP(I)-GAMd) 
CF  ( I )  =  -ZAZA/CGd) 


HP  =  -ZAZA  *  tGN  *  (  E*FKAP(I  )+GAM(I)  )  /  CGd) 

HP  =  HP  /  (  EFN  *  tGN  *  FKAPd)  ♦  ZAZA  ) 

HM  =  2.0  -  HP 
HM  =  SURT  (  HM  ) 

HP  S  bQRT  (  HP  ) 

Ml)  =  ZK  *  bQt  *  (  GAM  (I)  *  HM  -F  GNU  *  HP  ) 

Gd)  =  -ZA*F  (I)/CG(I) 

2/0  DGd)  =  Fd)*((1.0-CF(1)  )*tGN-CF(I)*EFN)/TUGAM(I) 
UFd)  =  -Gd  )*(  (1.0-CGd)  )*EFN-CG(I)*EGN)/TUGAM(I) 
300  CONI1NUE 
GO  10  309 
302  HAZ=bQRT  (ZA) 

I UAZ  =  2 •  0*ZA 
1ENU  =  1 
UO  30b  1-1 » K2M 

GAM ( 1 )  =  SQHr  (FKAP(1)**2-ZAZA) 

TUGAMd)  =  2.0  *  GAM  ( 1 )  ♦  1.0 
Mil  =  HAZ*(  rUAZ**GAM(l)  ) 

IF  (FKAPd/ .GT. 0.0)  GO  TO  304 

F (1)  =  -F  1 1 ) 

CF  (1 )  =  FKAPd)-GAM(l) 


CG(I) 


RAD00130 
RAD00140 
RAD00150 
RADOOlbO 
RAD00170 
RAD00180 
RAD00190 
RA000200 
RAD00210 
RAD00220 
RAD00230 
RAD00240 
RAOC0250 
RAD00260 
RAD00270 
RAD00280 
RAD00290 
RAD00300 
RAD00310 
RAD00320 
RAD00330 
RAD00340 
RAD00350 
RAD00360 
RAD00370 
RAD003b0 
RAD00390 
RAD00400 
RAD00410 
RAD0C420 
RAD00430 
RAD00440 
RAD00450 
RAD00460 
RAD00470 
RAD004H0 
RAD00490 
RAD00500 
RAD00510 
RAD00520 
RAD00530 
RAD00540 
RAD00550 
RAD005b0 
RAD00570 
RAD00580 
RAD00590 
RAD00600 
RADOOol 0 
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C6U)  =  -ZAZA/CM  I) 

AZKAP  =  CM1)/ZA 
GO  10  30b 

3UH  C6(l»  =  -FKAP(l)-GAMI) 

CMD  =  -ZAZA/CGUI 
AZKAP  -  -ZA/CGU) 

3US  G(I)  =  AZKAP*F II) 

OF  ID  =  “F  ( I )  * T U  AZ  / 1  OGAM  ( I ) 

UGU)  =  Z.U*F(I)*(1.0-CF(I)  J/TUGAHCI) 
3Ub  CONllNUt 

aus  ntab  =  1 

K1  -  LM-U 
M2  -  LM+2 

1UP  =  MAXU  (NK»NKP) 

00  aib  NO-1 # I UP 
Fa (NO 3  =  U.O 
FKPtNO)  =  U.U 
UFMHOI  =  U.U 
UFAPtNO)  =  O.U 
bFAlNU)  =  U.U 
bib  bFXP(NO)  =  U.U 
00  aib  N0=1#A2M 
bhtNUI  =  O.U 
bib  bG(NO)  =  O.U 
FI=U.UU  JB12b 
N  =  loU/H  ♦  .1 
KONt=l .0 

uo  aau  i-l fN 

abU  CALL  RAU1 

00  abb  1=1#K2M 

CALL  LOGoAM  ( TOGAM ( I ) *  U  # DUM 1  * 0UM2 ) 
RNORMU)  =  LXP  (-0UM1) 

Ml)  =  F  (I )  *RNOKM  ( 1 ) 

G(l)  =  G(I)*RNORMU) 

UF  (1 )  =  UFU)*RNORMtI) 

UGU)  =  UG  ( 1 )  *RNORM  ( 1 ) 
bF(l)  =  bF(I)*RNOHM(I) 
bGU)  =  bG(l)*HNORMU) 

CF 1 1 )  =  PKAPU) 
abb  CGU)  =  -FAAP(I) 

00  bUU  N=1*NA 
1  =  AGIN) 

FA(N)  =  FA(N)*RNORM(I) 

UFAIN)  =  UFA (N) *RNORM( I ) 
bUU  bFAIN)  =  bFA(N)*RNORMU) 

UO  blU  N-l * NKP 
1  =  AF  tN) 

FAPIN)  =  PAP(N)*RNORMli) 

UFAFIN)  =  UFAP(N)*RNOKM( I) 
blU  bFAPIN)  =  bFAPlN)*RNORMU  ) 

MRIlt  (  b»  1UU  ) 

MRlIt  lb»lUl)  Fl.RONt 
H  =  U . 12b 

XC'JI  =  AM1N1  IRCUT  #bb.U ) 

MR 1 1 L  l  b »  1U1  i  H»  XCUT 
NtM  =  1 
CALL  UtRlV 

NUON  =  IXCUf-l .U)/H+U.l 
UO  a/U  NUO-1 » NUON 
a/U  CALL  RKUT 
RX  =  XCUI 

IF  IRCUT .Lt.XCUI )  GO  TO  i78 
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RA000620 
RAD00630 
RAD00640 
RAD006S0 
RAD00660 
RA000670 
RA000680 
RA000690 
RAD00700 
RAD00710 
RAD00720 
RAD00730 
RAD00740 
RAD007b0 
RAD00760 
PAD00770 
RAD00780 
RAD00790 
RAD00800 
RAD00810 
RAD00620 
RA000830 
RA000840 
RAD00850 
RAD00860 
RAD00870 
RAD008B0 
RAD00890 
RAD00900 
RAD00910 
RAD00920 
RAD00930 
RAD00940 
RAD00950 
RAD009b0 
RAD00970 
RADOO90O 
RAD00990 
RAD01000 
RAD01010 
RAD01020 
RAD01030 
RAD01040 
RAD01050 
RAD01060 
RAD01070 
RADO1O0O 
RAD01090 
RAD01100 
RADO 1110 
RAD01 120 
RAD01130 
PAD01 140 
RAD01 150 
RADOllbO 
RAD01 170 
RADOllbO 
RAD01 190 
RAD01200 
RAD01210 
RAD01220 
RAD01230 


Nt*=l 

CALL  XDtKlV 
■  AVt  =  AMAA1  (KK#Q) 
lH-PI/MAVt 
H=FLOAT  (lH)/a.O 
IF  (H.GT.1.0)  H— 1  .  U 

OKI  It  (  6*  101  )  H»  KCUT 
LUb  =  (KCUI-KXI/HHO.l 
00  3/b  1-1 c LUb 
3/b  CALL  XHKUT 
3/B  CALL  WNOKM  ( KNURM t HX I 

IF  (  It NO  .to.  1  )  00  TO  390 

Ml  =  KM  ♦  £ 

M2  —  KM  ♦  3 
T  =  KK  *  HCUT 
CALL  bPHOtb(Y) 

00  38U  1  —  1 » K2M 
J1  =  LMK ( 1 )  ♦  1 
J2  =  LM1)  *1 

bltl)  =  NKt  *  bN(I)  *  fell)  *  d(Jl)  -  F ( I )  *  b(J2) 

CH  ( 1 )  =  HKt  *  G(I)  *  BIJ2)  ♦  SNd)  *  F  ( I  >  *  b(Jl> 

KNM  -  bOHT  (  bl(I)  *  bid)  ♦  CR(I)  ♦  CR ( 1 )  ) 
bid)  =  bid)  /  RUM 

CKd)  =  CHd)  /  RNM 

38U  CONI lNUt 
3SU  00  HUU  K-l t NK 
I  -  Kb  IK) 

L  TOl  =  LMKd*LKd)+LbtS(K)-l 

IF  l MOU (LTOI » H ) . Nt . 0 )  FK(K)=-FK(K) 

HUU  F K ( K )  =  F"  K  l  K  ) *HN0RM ( 1 ) 

00  HU 1  K-l #NKP 
1  =  KMK) 

L  T01  =  LMKd+LK ( 1 ) *LBb ( K ) —1 
IF  MUU(LrOT.H) .Nt.O)  FKP(K)=-FKP(K> 

HU  1  FKHIK)  =  FKP(K)*RNOKMI) 

HCU  CALL  HUM 
KETUKN 
tNU 


RAD01240 

RA001250 

RA001260 

RAD01270 

RAD01280 

RAD01290 

RAD01300 

RAD0131C 

RAD01320 

RAD01330 

RAD01340 

RAD01350 

RAD01360 

RAD0137C 

RAD01380 

RAD01390 

RAD01400 

RAD01410 

RAD01420 

RAD01430 

RAD01440 

RAD01450 

RA001460 

RAD01470 

RAD01480 

RAD01490 

RAD01500 

RA001510 

RAD01520 

RAD01530 

RAD01540 

RAD01570 

RAD01550 

RAD01560 

RA001580 

RAD01590 

RAD016U0 

RAD01610 


-  90  - 


SUBROUTINE  RKUT 


Purpose:  Performs  the  Runge-Kutta  integration.  The  routine  uses  indexed 

tabulated  values  of  the  bound-state  wavefunction  and  the  potential 
obtained  previously  by  interpolation  from  the  Waber  output. 

Radial  integration  is  performed  up  to  c  maximum  of  one-half  Bohr 
radius . 

Method:  Runge-Kutta  Integration  (Gill  Form).  A  fourth  order  integration 

scheme  in  which  the  error  in  each  step  is  of  the  order  h',  where 
h  is  the  interval  size. 


Subroutine  called:  DERIV 
Subroutine  called  by:  RADINT 

Labelled  Common:  DFUNC,  KUT,  LIMIT,  MAT,  TAPES 


Local 

Variables : 

Name 

Dimens  ion 

Mode 

Meaning 

J 

I 

Loop  index  on  Runge-Kutta  integration 

I 

I 

Loop  index  on;  number  of  free  wavefunctions 
number  of  matrix  elements  K  (hh  ' )  and 

K  (k'h)  1 

X 

Z 

R 

Incremental  variable  for  "small"  component 
of  free  electron  wavefunction 

ZP 

R 

Incremental  variable  for  "large"  component 
of  free  electron  wavefunction 

-  rl  - 


Return 


SlttMC  HKUto 

bUBKOUliNt  KKUT  RKU00010 

C  KUNGE-KUTTA  INTEGRATION  RKU00020 

COMMON/OFUNC/F (JO) *6(301 *DF (30) *06(30) *DFK (200) *DFKP<200> #CF(30) »HRKU00030 
COMMON  /KUI /RKl CHI »KK2 (4) »RK3 (4 ) »RK4(4) »K4 (4)  RKU00040 

C0MM0N/L1MH/JM»LM»KM»K2M»IEN0»NE<(»NK»NKP»JKB»LMKB»NTAB  RKU00050 

COMMON/MAT/bF ( 5P >  »bb<30) *FK (200) *FKP(200) •SFKI200) »SFKP(200) *RCUT  RKU00060 
COMMON  / I API  ./At] jUO) »bCF ( 1500) *FB( 1500) »G8( 1500) »GAMB» SCREEN  RKU00070 

00  150  J=1 1  ■*  RKU00060 

CALL  0ER1V  RKU00090 

00  120  1  =  1  ’<<.i  i  RKU00100 

Z  =  KM(J»*lOFtl»-KF.2(J)*SF(I)»  RKU00110 

ZP  =  HK1(J)*(UG(I) “NK 2 ( J ) *S6  d ) )  RKU00120 

Ml)  =  M1IWZ  RKU00130 

Gtl)  =  Gd)*H*2P  RKU00140 

bFd)  =  bF  d ) *3*  0*Z-KK3( J) *DF ( 1 )  RKUOOlbO 

1ZU  SGd)  =  bG(l)+3.0*ZH-HK3(J)«DG(I)  RKU00160 

00  1Z1  1=1 »NK  RKUOC170 

Z  =  NK11J>*(0FM1)-KK2(J>*SFK(1M  RKU00180 

FK(1)  =  FM1>*FI*Z  RKU00190 

121  bFKIl)  =  bFK  d ) ♦3«0*Z“HK3( J) *DFK ( I )  RKU002U0 

00  1Z2  1=1 »NKP  RKU00210 

Z  =  KM(J)*tOFKPtl)-KK2(J)*SFKPU)  )  RKU00220 

FKPd)  =  FKP ( 1 ) ♦H*Z  RKU00230 

1ZZ  bFKPd)  =  bFKPd ) ♦3«0*Z“KK3( J)*DFKPl I )  RKU00240 

NfAb=Nl Ab+KH ( J)  RKU00250 

130  CONtlNUE  RKU00260 

tfEIUNN  RKU00270 

tNO  RKU00280 


i 
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SUBROUTINE  S1NDEX 


Purpose:  To  catalogue  and  index  matrix  elements  consistent  with  selection 

rules . 

Method:  Routine  examines  trie  angular  momentum,  values  to  find  and  index  all 

the  matrix  elements  compatible  with  the  selection  rules  and  with 
the  input  cutoff  values  of  the  quantum  numbers  and  records  the 
free  electron  h  and  the  photon  l  values  for  each  matrix  element 
in  reference  vectors  (separately  for  the  two  kinds  of  matrix 
elements ) . 

Subroutine  cal le d :  None 
Subroutine  called  by:  PELEC 
Labelled  •*  ;;  :  QUANT,  VECT ,  LIMIT 

Local  /ariables : 


Dimension 


Me an in 


Loop  index  on  max  k  for  electron 


Loop  index  on  number  of  free  electron 
states 


LSUM  I  l  +  l  , 

K  -K 

LD1FF  I  I  4  -l  ,  I 

|  H  -H  I 

L  I  Photon  angular  momentum 
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Selection  rule  violated 


SiBHC  SINU 

SUBHOU! INt  blNDLX 

COMMON/ QUAN  f /LK ( 30 ) *LMK(3Q) # JK<30> »FKAP(30> »SN(30) *SI (30) »CR(30> 

COMMON/ VECT/KF (200) »KG(2U0) »LbES(200) »LBS(20Q) *LKB 

COMMON  /L1M1T/JM»LK»KM»K2M#IEND#NEW#NK»NKP»JK8*LMKB»NTAB 

NK'O 

NKF-U 

DO  b  1=1*KM 
*=2*1-1 
J-K  ♦  1 

FKAP(K)=-1 
FKAFl J)=l 
LK (K ) =1-1 
LMU)=I 
LMK l K ) - 1 
LMM  J)=l-1 
SN(K)=-1.U 
bNl  J)=1.U 

JK(K)=LMK)+LMMK) 
b  JK(J)=JK(K) 

00  1U  K=lrK2M 
LbUM  =  LK(K)+LMKO 
IF  l MOL) (LbUM #  2 )  )  3b*3U»3b 
3U  L-U 

00  10  40 
Ob  L-l 

bb  IF  IL.OT.LbUM)  00  TO  4b 
4U  LU1FF  =  I  Abb  (LK(K)-LMKb) 

IF  (L.LT.LUlFF)  00  TO  45 

NK-NK+1 

LbtblNK)=L+l 

KG(NK)=K 

4b  LT=LMMK)+LKb 

IF  IL.OT.LT )  00  TO  bU 

LU=lAbb(LMMK)-LKb) 

IF  IL.L1 .LU)  GO  TO  bO 

NKF— NKF+1 
LbblNKK)=L+l 
KF (NKF)=K 
bU  L=L+2 

IF  IL.Lt.LM)  00  TO  bb 

10  COfJUNUt 
KLTUKf) 
tNU 


SND00010 

SND00020 

SND00030 

SN000040 

SND00050 

SND00060 

SN000070 

SN000080 

SND00090 

SND00100 

SND30110 

SN000120 

SND00130 

SND00140 

5N000150 

SND00160 

SND00170 

SNDOOlbO 

SND00190 

SND00200 

SND00210 

SND00220 

SND00230 

SND00240 

SND00250 

SND00260 

SND00270 

SN000260 

SND00290 

SND00300 

SND00310 

SND00320 

SND00330 

SND00340 

SND00350 

SND00360 

SND00370 

SND00380 

SND00390 

SN000400 

SND00410 

SND00420 

SND00430 


SUBROUTINE  3 FREES 


Purpose:  Computes  the  values  of  the  spherical  Bessel  function. 

Method:  The  zero  order  function  is  obtained  as  sir.  R/R.  For  small 

argument  the  higher  order  functions  are  calculated  from  the 
power  series  expansion.  For  large  argument,  the  first  order 
function  is  computed  from  its  explicit  sinusoidal  representation 
and  the  remaining  ones  obtained  by  recursion  relations.  For 
intermediate  argumei . *•  ~ ,  the  lower  order  functions  are  obtained  by 
recursion,  the  higher  order  by  nower  series. 

Subroutines  called:  None 

Subroutine  called  by:  DERIV,  XDERIV,  RADINT 

Labelled  Common:  BESSEL 

Argument  sequence:  (R) 

Argument  List: 

Name _ Dimens  ion _ Mode _ Meaning _ _ 

R  R  Bessel  function  argument;  photon  momentum 

*  radius  from  DERIV,  XDERIV;  free  electron 
momentum  x  radius  from  RADINT. 


Local  Variables : 

Name _ Dime  ns  ion _ Mode 

NR  I 

L  I 

SER  R 

TER  R 


Meaning _ 

Defines  transition  from  power  series  to 
recursion. 

Loop  index  on  calculating  Bessel  function. 
Partial  sum  of  power  series  expansion. 

Last  term  of  series  expansion. 


-  QQ  _ 


Calculate  Bessel  function 


*1BML  SPHB 

bUdKOUTlNt  SPHdtS  (H) 

L  bPHtKKAL  dtSbtt  FUNCTION 

COMMUN/btbbtL/FLdb)  #PC  ( lb)  »OF ( 15#  lb)  »M1  »M2»B(  15) 


d(] 

L) 

=  b 

IN  IK ) 

• 

OR 

IF 

( 

Ml 

.to. 

1 

) 

GO 

TO 

16 

NK 

- 

H  ♦ 

2.0 

IF 

1 

NK 

.LI  . 

2 

) 

GO 

TO 

12 

B(2> 

=  ( 

H  ( 1 ) 

- 

COb(K) 

)  •  OR 

IF 

l 

Ml 

•  tu. 

2 

) 

GO 

TO 

16 

IF 

( 

NK 

.61 . 

Ml 

) 

NK=M1 

UO 

11 

L 

=  2» 

NK 

11  B(L  +  1 )  =  FL(L)  *  B(L)  *  OR  -  b(L-l) 

IF  l  NK  -  Ml  I  13.  16*16 

12  NN  =  2 

lb  MAb  =  U.b  •  H  *  R 
UO  lb  L  -  NN»  M2 

j  =  u 
btK  =  1.0 
ltK  =  1.0 
1H  J  :  J  4  l 

ltK  =  -  TtK  *  HAS  *  OF(L»J) 
btK  =  btK  ♦  TtK 

IF  1  AbS(TtK)  .GT.  (  Q.OU0001  *  ABS(SER)  )  > 
lb  B(L)  =  SEK  *  PC (L)  *  (  K**(L-1)  ) 
lt>  KtlUNN 

tuu 


GO  TO  14 


SPH00010 

SPH00020 

SPH00030 

SPH00040 

SPH00050 

SPH00060 

SPH00070 

SPH00080 

SPH00090 

SPH00100 

SPH00110 

SPH00120 

SPH00130 

SPH00140 

SPH001S0 

SPH00160 

SPH00170 

SPHOOlttO 

SPH00190 

SPH00200 

SPH0C210 

SPH00220 

SPH00230 

SPH00240 

SPH00250 

SPH0G260 

SPH00270 
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SUBROUTINE  WN$RM 


Purpose:  Computes  the  normalization  factors  for  the  free-electron  wave- 

function  after  termination  of  the  numerical  integration. 

Method:  The  normalization  factors  are  determined  by  a  matching  to  the 

W.K.B.  solution.  Derivatives  of  the  potential  in  the  W.K.B. 
expression  are  obtained  from  a  numerical  fit  of  an  exponential 
to  the  screening  function  (ratio  of  screened  to  unscreened 
potential)  at  the  cutoff  radius.  When  the  W.K.B.  conditions 
fail  (vanishing  kinetic  energy  and  high  angular  momentum),  the 
normalization  factor  of  the  previous  wavefunctions  is  used. 
Subroutine  called:  None 
Subroutine  called  by:  RADINT 

Variables  in  unlabelled  Common:  PI,  HALFPI,  F^URPI,  RAD,  SQ2,  Q,  ZA,  ZAZA, 

EFN,  EGN,  V,  CG,  GAM 
Labelled  common:  DFUNC,  LIMIT,  TAPES 
Argument  sequence:  (HN  RM,  RX) 

Argument  List: 

Name _ Dimens  ion _ Mode _ Meaning _ 

HNs2>KM  30  R  Normalization  factors 

RX  P  Terminal  radius  on  numerical  integration 


Local  Variables: 

Name _ Dimension _ Mode 

0VERR  R 

EMVPl  R 

EMV  R 

0RIGH  R 

I  I 


Meaning _ _ 

Reciprocial  of  terminal  radius 

Free  electron  energy  -  potential  energy  +  1 

Free  electron  energy  -  potential  energy 

Radial  grid  interval 

Loop  index  on  number  of  free  electron 
states 
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RA,  RB 


EP 

P'Q 

P3QP 


R  Tabulated  values  of  screening  factor 

at  grid  points  bracketing  terminal  radius 

R  Screening  attenuation  coefficient 

R  ?  in  W.K.B.  solution 

R  Derivative  of  Pc 
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SittMC  WNUR 

bUBHOUTINt  WNOKM  <HNORM#RX)  MNR00010 
COMMON  P1»hALFPI#FOURP1#RAD*SQ2#Q#ZA#ZAZA#EFN.EGN»V.CG(30)*GAM(30)WNR00020 
COMMON/UFUNC/F<30) #0(30) #DF(30>  #DG«30) #DFK(200) #DFKP(200) »CF(30> »H*NR00030 


COMMON/  L I  M I T/  JM » EM  #  KM » K2F*  ♦  I  END  »  NEW  *  NK  #  NKP  #  JKB  »  LMKB  »  NT  AB 
COMMON/TAPtb/X(lbOU) #SCF (1S00) #FB(1500) »GB(1500) »SAMB»SCREEN 
UIMtNblON  NNOKM ( 30 ) 

OVEKK  s  1,0  /  RX 
OVKHbO  =  OVERK*OVEKH 
tMVPl  =  EGN-V 
LMV  =  EMVP1-1.U 
Y  =  -V/EMVP1 

ORlGH  =  X ( NT Add )  -  X(NTAB) 

00  3  1-1  #K2M 
PK  =  CF  (  1  ) 

FKbG  =  PK* ( PKd , 0 ) 

FKR  =  PR  *  OVEKH 
FOR  =  F  KbO  *  OVRRbO 
KA  =  bCF (Nl Ad> 

Kd  =  bCF ( N I Ad+ 1 ) 

EP  =  ALOG(KA/KO)  /ORIGH 
VI  =  EP  ♦  OVEKK 
V2  =  VI  *  VI  ♦  OVRRbO 
V3  =  VI  *  V2  ♦  2.0  *  OVRRbQ 
VlX  =  V  1*Y 
V2X  =  V2*Y 
VXbU  =  VlX  *  VlX 
FKVK  =  V 1  X*F  KK 
PbU  =  EMV  *  EMV  ~  1 
IF  (PbG.LE.OVKKbQ) 

PbUP  =  2.0  *  EMV  *  V 


*  (VI  ♦  OVERR) 


0  ♦  FKVR  -  0.75  *  VXSO  ♦  0.5  *  V2X  -  FOR 
GO  TO  1 

*  VI  ♦  (FKK  -  1.5*V1X)  *  VXSQ  +  2.0*  VlX* 


1  -  U.b  *  V3  *  Y  -  OVEKR  *  (FKVR  ♦  PK*V2X  -  2.0*F0R) 

IEKM3  =  0 ( 1 )  *  ( (J ■  b*VlX— FKR+U .25* (PSQP/PSQ)  )  +F  C 1 ) *EMVP1 
A  =  (Pl/(boKT  (PbQ) *EMVP1 ) ) * (PSQ*G( I )  *G(I>  +TERM3*TERM3) 
KNOKMU)  =  l.U/SQRT  (A) 

GO  IU  3 

IF  U.OT.l)  GO  TO  2 

KNQKM  ( 1 )  -bCKEtN 
GO  IU  b 

KNUKrfl 1 )=RNOHM< 1-1 ) 

CON  I  IhiUE 
keiukn 

ENO 


WNR00040 
MNR00050 
MNR00060 
WNR00070 
WNR000B0 
WNR00090 
WNR00100 
WNR00110 
WUR00120 
WNR00130 
WNR00140 
WNR00150 
WNR00160 
WNR00170 
WNROOlflC 
WNR00190 
WNRG02U0 
WNR0021 0 
WNR00220 
WNR00230 
WUR00240 
WNR00250 
WNR00260 
WNR00270 
WNROO2B0 
WNR00290 
V2XWNR00300 
WNR00310 
WNR00320 
WNR00330 
WNR00340 
WNR00350 
WNR00360 
WNR00370 
WNR003B0 
WNR00390 
WNR00400 
WNROO410 
WNR00420 
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SUBROUTINE  XDERIV 


Purpose:  Supplies  the  derivatives  for  the  Runge-Kutta  integration  of 

the  hound-state  wavefunctions  and  the  integrands  of  the  matrix 
elements  for  numerical  integration  beyond  one-half  Bohr  radius. 
Method:  Calculates  the  derivatives  of  the  radial  components  from  the 

coupled  Dirac  radial  equations  and  the  integrand  of  the  matrix 
elements  using  values  obtained  by  linear  interpolation  on  the 
Waber  grid  of  bound  state  wavefunctions. 

Subroutines  called:  SPHBES,  INTERP 

Subroutines  called  by:  XRKUT,  RADINT 

Labelled  Common:  BESSEL,  DFUNCT,  LIMIT,  0NWARD,  VECT 


Local  Variables: 

Name  Dimension 

Mode 

Meaning 

Z 

R 

Photon  momentum  *  radius 

N 

I 

Indexing  variable 
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S1BMC  XDtH 

SUBKOUlINt  XDtH IV 

C  COMPurtb  DERIVATIVES 

COMMON  PI tHALFPl tFOUKPI »RA0»SU2»Q»ZAt ZAZA»EFN»EGN»  V»CG( 30 
COMMON/BtSbtL/FL(lb)  »PC(lb)  »0F(15»15)  »M1»M2»BU5) 
COMMON/DFUNC/F  (30)  >G(30)  #DF(30>  #06(30)  »DFK(200)  *DFKP(200) 
COMMON/ L I M i I / JM » LM » KM » K2M *  I END > NL W » NK  * NKP » JKB * LMKB » NT  AB 
COMMON  /ON»ANO/RX#SCX»GBX»FBX 

COMMON/ VtC I /KF  (200)  #KG(200)  »LbES(200)  »LBS(200)  »LKB 

IF  (NtW.tQ.O)  GO  TO  30 

NEW  =  0 

CALL  INTLHP 

V  =  -  SCX  /  KX 

Z=U*KX 

CALL  SPHBtS  (Z) 

00  10  31 

30  NE*  =  1 

31  00  3b  N=1*K,2M 

UF (N)  =  CF(N)*F (N)/KX-<tFN-V) *G(N) 

3b  UG ( N )  =  CG(N)*G(N)/KX-MEGN-V)*F(N) 

00  HO  N- 1 » NK 
1  =  KG(N) 

L  =  LbtS(N) 

HO  UFK(N)  =  b(L)*G(  I )  -*FbX 
00  HI  N— 1 » NKP 
1  =  (N) 

L  =  LbS(N) 

HI  OFKP(N)  =  b(L) *F ( I ) *ObX 
HEIUKN 
tNO 


XOROOOlO 
XDR00020 
) » GAM (30) XDR00030 
XDROOOV  0 
»CF(30) »HXDROOOSO 
XDR00060 
XDR00070 
XDR00080 
XDR00090 
XDROOlOO 
XDR00110 
XDR00120 
XOR00130 
XDR00140 
XDR00150 
XDR00160 
XDR00170 
XDR00180 
XDR00190 
XDR00200 
XDR0021 0 
XDR00220 
XDR00230 
XDR00240 
XDR00250 
XDR00260 
XDR00270 
XDR00280 
XDR00290 
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This  page  is  intentionally  left  blank. 
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SUBROUTINE  XRKUT 


Purpose:  Perfonas  the  Runge-Kutta  integration  for  radial  values  beyond 

one-half  Bohr  radius.  The  routine  uses  linearly  interpolated 
values  of  the  bound-state  wave functions  and  the  potential. 
Method:  Runge-Kutta  Integration  (Gill  Form) 

Subroutine  called:  XBERIV 

Subroutine  calle  a  oy:  RADIHT 

Labelled  Oonmon:  DFUNC,  KUT,  LIMIT,  MAT  0NWARD 

Local  Variables : 


-  Ill  - 


This  routine  is  used  for  Runge- 
Kutta  integration  beyond  one-half 
Bohr  radius 


Call  XDERIV  tc 
calculate  derivatives 


Loop  for  Runge- 
Kutta  integration 
to  £.30 
V. ' 

Loop  on  wave 
functions  to 
020* 


Loop  on  matrix 
elements 

K  (ax’)  to/l21 

l  V  y> 


Loop  on  matrix 
elements 

K  (h’k)  to  122 

«  ./ 


Return 
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SIBMt  XKAU 

SOBKOOlINt  XKKUT  XRK00010 

C  KUNGE-KUT1A  INTEGRATION  XRK00020 

LOMMON/UFUNC/F  (30)  »G(30)  »DF(30)  »DG(30)  »DFK(200>  »DFKP(200)  »CF(30>  »HXRK00030 
COMMON  /RUf/  KM(4)  »KK2(4)  »R*3(4)  »RK4(4)  ,K4{4)  XRK00040 

COMMON/LlMir/JM»LM»KM»K2M»IENU»NEW»NK»NKP»JKB»LMKB»NTAB  XRK00050 

COMMON/MAf/SF (30) »SG( 30 >  »FK (200 ) »FKP(200) »SFK(200) »SFKP(200) »RCUT  XRK00060 
COMMON  /ONWARD/ KX » SCX » GtiX » FBX  XRK00070 

00  130  J=l»4  XRKOOOBO 

CALL  XUEK1V  XRK00090 

00  120  1=1»K2M  XRK00100 

L  =  HKl(J)*(OF(I )“HK2( J) *SF ( I ) )  XRK001 10 

ZP  =  KKl(J)*(UG( I)-KK2( J)*S6(I> )  XRK00120 

Ml)  =  Mll*H*Z  XRK00130 

G ( 1 )  =  G( 1 > *H*ZP  XRK00140 

SFU)  =  SF  ( 1 )  ♦3«0*Z“HK3( J) *DF  ( I )  XRK00150 

120  SGI  1 )  =  SG(l>+3.0*/H-RK3(J)*DG<n  XRKOOlbO 

00  121  I  —  1 » NK  XRK00170 

Z  =  KR1(J)*«0FMI)-HK2(J)*SFK(I))  XRKOO10O 

FK i 1 )  =  FK(X)+M*Z  XRK00190 

121  SFK(l)  -  SFK ( 1 ) ♦3»U*Z“KK3t J) *0FK ( I )  XRK00200 

OO  122  1=1»NKP  XRK00210 

Z  =  KK1 ( J) * (OFKP ( I ) -KK2 l J) *SFKP< 1 ) )  XRK00220 

FKPll)  =  FKHII)+H*Z  XRK00230 

122  SFKP(l)  =  SFKP( 1 ) ♦3»U*Z“HK3( J) *DFKP ( I )  XRK00240 

KX  =  KX*KK4(J)*H  XRK00250 

10U  CONIlNUt  XRK00260 

KtIUKN  XRK00270 

tNU  XRK00280 
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